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Preface 


The  product  of  this  thesis  (MMSOFE)  should  prove  to  be  a  useful  software  tool 
to  the  US  Air  Force  guidance  and  control  community  and  a  valuable  starting  point  for 
developing  other  multiple  model  application  software.  MMSOFE  could  prove  to  be  either 
a  boon  or  a  bain  to  future  AFIT  guidance  and  control  students  who  may  be  required  to  use 
MMSOFE  to  perform  multiple  model  adaptive  estimation  with  N  elemental  Kalman  filters 
instead  of  just  a  simple  single  Kalman  filter  simulation  like  students  in  previous  years.  To 
those  future  students,  I  say,  “You’re  welcome.”  While  originally  intended  as  a  prelimary 
software  development  portion  of  this  thesis  to  be  used  on  a  larger  more  complex  problem, 
the  development  of  MMSOFE  itself  proved  to  be  easily  thesis-worthy  and  challenging. 

First  and  foremost,  neither  the  successfully  completion  of  this  thesis  nor  my  enduring 
the  entire  AFIT  experience  would  have  been  possible  without  the  excellent  support  of  my 
loving  wife,  Vera,  and  my  dear  children;  Lynsa,  Eric,  Amber,  and  Lauren  (15,  12,  9,  and 
5  years).  Their  unselfish  sacrifice  of  time  with  me  and  of  my  attention  speaks  magnitudes 
for  their  characters  and  for  their  continuing  loving  and  prayerful  support.  Such  cannot  be 
repaid,  nor  should  the  sacrifice  of  any  AFIT  student’s  family  be  taken  for  granted.  To  my 
family,  “I  love  you,  and  thank  you  for  enduring  with  me.” 

My  gratitude  especially  goes  to  Dr.  Peter  May  beck  for  his  painstaking  proofreading, 
for  his  advice  and  direction,  for  being  a  touchstone  for  sanity-checking,  and  most  especially 
for  his  continuing  prayers  of  support.  Additionally,  LtCol  Rx>bert  Riggins  deserves  a  special 
“thank  you”  for  discussing  software  bugs  and  development  problems  and  for  his  excellent 
comments  and  corrections  to  this  thesis.  My  appreciation  is  also  due  Maj  Randall  Paschall 
for  taking  time  in  his  overfull  schedule  to  proofread  and  to  those  at  Wright  Laboratory 
for  bearing  with  me  while  I  finished  this  effort.  Last,  but  not  least,  “Thank  you”  to  Mr 
Donald  E.  Smith  of  Co-Operative  Engineering  who  spent  many  hours  ensuring  adequate 
disk  space  was  available  and  for  keeping  “MY”  SPARC  operating. 

p.s.  Never  again! 


Robert  L.  Nielsen 
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Abstract 

Multiple  model  Kalman  Filter  (KF)  techniques  are  used  extensively  for  Multiple 
Model  Adaptive  Estimation  (MMAE),  Multiple  Model  Adaptive  Control  (MMAC),  and 
Distributed  Kalman  Filter  (DKF)  applications  to  determine  Bayesian-blended  optimal 
estimates  of  states,  uncertain  parameters,  and  optimal  control  signals.  Multiple  model 
methods  are  used  for  sensor  management.  Failure  Detection  and  Isolation  (FDI),  and  other 
Guidance  and  Control  (G&C)  applications.  A  simulation  tool  called  the  Multiple  Model 
Simulation  for  Optimal  Filter  Evaluation  (MMSOFE)  has  been  in  this  research.  MMSOFE 
is  based  on  the  well-benchmarked  single  Kalman  filter  tool  called  Multimode  Simulation  for 
Optimal  Filter  Evaluation  (MSOFE).  MMSOFE  is  a  highly  portable  and  versatile  multiple 
and  single  Kalman  filter  evaluation  tool.  It  is  capable  of  performing  simulations  with  one 
filter  or  up  to  98  elemental  filters  in  a  multiple  model  adaptive  filter  structure.  It  can 
be  adapted  easily  for  other  multiple  model  applications.  MMSOFE  was  applied  to  failure 
detection  and  isolation  of  measurement  jamming-  and  spoofing-type  failures,  similar  to 
jamming  and  spoofing  of  a  Global  Positioning  System  (GPS).  A  satellite  orbit  estimation 
test  case  was  used. 


DEVELOPMENT  OF  A  PERFORMANCE  EVALUATION 
TOOL  (MMSOFE)  FOR  DETECTION  OF  FAILURES 
WITH  MULTIPLE  MODEL  ADAPTIVE  ESTIMATION  (MMAE) 

/.  Introduction 

Military  and  civilian  aircraft  employ  a  wide  variety  of  Inertial  Navigation  Systems 
(INS),  Global  Positioning  System  (GPS)  receivers,  and  numerous  other  navigational  aids. 
While  each  system  can  function  independently,  the  greatest  advantage  is  gained  through 
integrating  the  systems  available  on  a  particular  aircraft.  This  integration  is  frequently 
accomplished  using  a  Kalman  Filter  (KF).  With  improved  position  and  velocity  estimates 
obtained  with  the  KF-based  integrated  systems,  both  the  purview  of  applications  and 
system  precision  increase  substantially.  A  few  of  these  integrated  system  applications 
are:  precise  navigation  of  preplanned  flight  trajectories,  truly  automated  landing  approach 
systems,  air  refueling,  formation  flying,  possible  aiding  of  carrier  landing,  and  superior 
weapon  delivery  (be  it  air-to-air,  air-to-surface,  surface-to-air,  or  surface-to-surface). 

The  design  and  analysis  of  KF-based  integrated  system  designs  has  mostly  been 
accomplished  using  problem-specific  computer  simulation  software.  While  this  problem- 
specific  simulation/analysis  method  is  functional,  generating  the  model-independent  KF- 
common  code  is  a  time  consuming,  frequently  redundant  effort  and  further  enhances  the 
opportunity  for  error  -  maybe  even  "hidden”  errors  which  are  not  readily  observable  and 
could  go  undetected  during  development.  Alternatively,  analysis  of  system  designs  involv¬ 
ing  either  linear  or  extended  KFs  has  frequently  been  accomplished  by  running  computer 
simulations  using  a  program  called  Multimode  Simiilation  for  Optimal  Filter  Evaluation 
(MSOFE)  (24).  MSOFE  provides  a  software  tool  for  analysis  of  KF  designs  involving 
single  KFs.  It  allows  for  Monte  Carlo  simulations  of  linear  or  extended  Kalman  Filters, 
covariance  analysis  of  linear  or  linearized  filters,  or  both  types  of  analysis  simultaneously. 
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A  primary  import  of  MSOFE  is  that  a  complete  KF  analysis  is  possible  with  MSOFE  by 
using  its  Monte  Carlo  run  simulation  capability  without  conducting  flight  tests. 

However,  although  MSOFE  is  an  excellent  tool  for  simulating  and  analyzing  designs 
of  single  KFs  without  “reinventing-t he- wheel”  for  each  problem,  it  was  not  intended  to 
simulate  systems  which  utilize  multiple  KFs,  as  required  in  Multiple  Model  Adaptive  Es¬ 
timation  (MMAE),  Multiple  Model  Adaptive  Control  (MMAC),  or  Distributed  Kalman 
Filter  (DKF)  techniques.  MMAE,  MMAC,  and  DKF  are  discussed  in  more  detail  in  Sec¬ 
tion  1.5.2,  and  the  motivation  for  developing  a  Multiple  Model  Simulation  for  Optimal 
Filter  Evaluation  (MMSOFE)  is  discussed  in  Sections  1.1  and  1.2. 

1.1  Background 

MMAE  has  been  used  for  many  applications,  including:  Failure  Detection  and  Iso¬ 
lation  (FDI)  of  aircraft  control  surface  failures  (8)  and  reconfiguring  the  flight  control 
system  to  counter  the  failure  effects  in  such  situations  (33),  variable  stability  aircraft  flight 
control  (21),  forward-looking  airborne  target  tracking  system  (26),  MMAE  and  moving- 
bank  MMAE  for  sensing  and  controlling  vibrations  in  a  flexible  space  structure(7,  11), 
and  adaptive  control  of  a  robot  arm  (29),  to  name  a  few.  Another  prime  application  for 
MMAE  relates  to  research  being  conducted  by  the  Central  Inertial  Guidance  Test  Facility 
(CIGTF),  6585th  Test  Group,  Air  Force  Systems  Command  (AFSC),  Holloman  AFB,  NM 
to  determine  the  vulnerability  of  the  GPS  to  jamming  and  spoofing  (36). 

This  latter  application  was  the  original  impetus  for  this  thesis.  MMAE  can  be 
utilized  to  detect  and  isolate  jamming  or  spoofing  of  a  GPS  or  the  total  loss  of  a  GPS 
pseudorange  input  to  the  receiver.  Jamming  can  be  simplistically  described  as  bombarding 
the  GPS  receiver  with  broad-band  electromagnetic  noise.  In  the  case  of  jamming,  there  is 
no  eflfort  made  to  mislead  the  GPS  system  into  calculating  erroneous  results,  but  rather 
the  intent  is  to  prevent  reception  of  the  GPS  satellite  signals  entirely.  In  contradistinction, 
spoofing  is  employed  with  the  intentional  purpose  of  providing  misleading  information  to 
the  GPS  system  and  thereby  causing  the  system  to  calculate  erroneous  results.  Since 
effective  spoofing  must  be  clandestine,  the  spooler  attempts  to  mimic  the  GPS  satellite 
signal  while  adding  slight,  misleading  modifications  to  the  original  signal.  Undetected  and 
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Figure  1.1.  Types  of  Figures 

uncountered,  this  would  cause  a  GPS  -dependent-user  to  be  drawn  away  from  his  intended 
destination  due  to  the  resulting  navigation-related  errors.  Lastly,  another  type  of  failure  is 
merely  the  total  loss  of  GPS  pseudorange  measurement  data  without  any  hostile  intent,  but 
possibly  due  just  to  simple  antenna  shielding  or  signal  attenuation.  Serious  degradation 
in  performance  in  the  applications  mentioned  previously  would  be  the  final  result.  While 
the  errors  caused  by  all  three  types  of  failure  may  appear  very  similar,  it  is  important  to 
identify  which  of  these  three  types  occurred  and  caused  the  errors.  These  types  of  errors 
are  depict'-^  in  Figure  1.1  and  generally  represent  the  types  of  failures  simulated  to  occur 
in  the  orbit  range  measurements  for  the  orbit  estimation  problem  (to  be  described  in  detail 
in  Section  3.2)  used  as  the  test  case  for  MMSOFE  development. 

Figure  1.1  depicts  two  types  of  failures  without  specifying  the  type  of  desired  signal. 
With  regard  to  this  figure,  the  desired  signal  could  represent  a  physical  quantity  such  as 
distance,  velocity,  acceleration,  or  impulse  or  even  a  value  such  as  noise  variance.  For  the 
orbit  estimation  problem,  the  desired  signal  represents  either  the  orbit  range  measurement 
or  the  orbit  range  measurement  noise  variance.  After  a  specified  time,  the  range  measure- 
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ment  signaJs  are  corrupted  by  adding  a  bias  for  the  jump  failure  and  adding  a  ramped 
offset  for  the  ramp  failure.  In  the  case  of  the  measurement  noise  variance,  the  variance  is 
adjusted  by  adding  a  bias  or  ramp  at  the  time  of  failure. 

In  the  GPS  problem,  the  same  methods  could  be  used  for  the  failure  injection  (see 
Section  1.1).  The  bias  and  ramp  in  the  measurement  noise  variance  would  correspond  to 
the  jamming.  However,  while  a  jump  in  the  measurement  itself  could  relate  to  an  entire 
loss  in  the  pseudorange  signal,  either  a  jump  or  ramp  in  the  measurement  might  be  used 
to  simulate  simple  or  intelligent  spoofing,  respectively. 

Additionally,  the  failure  thresholds  shown  in  Figure  1.1  are  not  readily  known  a  priori. 
The  thresholds  could  be  set  a  priori  and  intelligently  if  a  great  deal  were  known  about  the 
trajectory  prior  to  flight  and  if  false  failure  alarms  and  missed  alarms  could  be  tolerated. 
A  more  realistic  and  practical  method  would  be  to  use  the  real-time  filter-computed 
covariances  of  the  residuals  to  determine  the  scalar  standard  deviations,  or  l-sigma  values 
(square  root  of  the  variances).  Some  multiple  of  these  l-sigma  values  could  then  be  used  as 
threshold  values  for  the  residuals  to  detect  failures  through  residual  monitoring  or  multiple 
model  hypothesis  test  methods  (to  be  discussed  later  in  more  detail). 

1.2  Problem  Definition 

The  problem  is  to  develop  a  simulation  and  analysis  software  tool  to  aid  controls 
designers  and  system  integrators  and  enable  them  to  concentrate  on  problems  specific  to 
multiple  model  estimation  algorithms,  by  eliminating  the  “reinvention”  of  analysis  tool 
software  and  thus  minimizing  the  sources  of  errors  resulting  from  the  “reinvention”  effort 
of  the  software  development  for  eacn  design  application.  Since  MSOFE  (4),  (24)  is  a  widely 
known,  benchmarked  program,  it  was  chosen  as  a  basis  upon  which  to  build,  rather  than 
building  a  new  independent  stand-alone  MMSOFE.  This  will  make  it  much  easier  for  an 
MSOFE-familiar-user  to  run  MMSOFE  without  any  great  learning  curve. 

Therefore,  the  primary  goal  of  this  thesis  is  to  modify  MSOFE  to  run  MMAE  algo¬ 
rithms.  This  modified  program  will  be  referred  to  as  Multiple  Model  Sim’-  ation  for  Op¬ 
timal  Filter  Evaluation,  or  MMSOFE.  Due  to  the  adaptive  nature  of  MMAE  algorithms. 
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only  Monte  Carlo  performance  evaluations  are  pursued,  and  there  are  no  analytically  cor¬ 
rect  equations  for  the  covariance  analysis  of  an  MMAE.  The  Satellite  Orbit  Estimation 
problem  described  in  detail  by  Maybeck  (19:pp.  46-48)  and  in  the  MSOFE  users’  manual 
(24)  will  be  used  for  the  test  problem  during  MMSOFE  development.  To  this  basic  prob¬ 
lem  will  be  added  the  simulated  failures  discussed  in  the  previous  section  to  demonstrate 
MMAE  performance  attributes.  MMSOFE  will  have  application  to  a  wide  variety  of  prob¬ 
lems  in  the  Air  Force  and  in  industry.  MMSOFE  will  initially  be  developed  to  produce 
optimal  blended  state  estimates  and  covariances,  to  calculate  parameter  estimates,  and  for 
failure  detection  and  isolation  (FDI).  Follow-on  goals  are  to  modify  MMSOFE  further  in 
order  to  permit  MMAE-based  control  or  MMAC,  to  investigate  adapting  MMSOFE  for 
systems  requiring  Distributed  Kalman  Filters  (DKF)  as  well,  and  to  apply  MMSOFE  to 
a  larger,  more  complex  problem  such  as  FDI  in  a  GPS/INS  integrated  system. 

1.3  Scope 

The  failure  types  considered  in  this  research  will  be  simulated  either  as  a  jump  or 
ramp  in  either  the  system  measurement  value  or  the  system  measurement  noise  variance.  A 
jump  fmlure  in  this  context  consists  of  a  nearly  instantaneous  change  in  the  signal  offset  or 
noise  variance  while  a  ramp  failure  is  a  gradually  increasing  change  in  the  bias  measurement 
offset  or  noise  variance.  Unlike  many  FDI  algorithms,  it  is  unnecessary  to  know  and  model 
the  exact  magnitude  of  these  jumps  in  advance,  since  the  MMAE  algorithm  can  estimate 
their  magnitude  automatically  if  various  hypothesized  jump  values  (hypothesized  ramps 
will  not  be  required)  are  modeled  in  the  elemental  filters.  Further,  while  uiese  type  failures 
may  be  applicable  to  many  different  systems,  the  primary  purpose  of  this  research  is  to 
develop  the  MMSOFE  tool  itself.  The  criteria  for  completion  of  the  research  are  presented 
in  Section  1.6. 

1.4  Assumptions 

The  assumptions  given  below  provide  the  scenario  for  the  MMSOFE  development 
and  the  simulation  used  in  this  research.  The  reader  can  further  investigate  their  import 
in  the  referenced  documents. 
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1.  Computer  simulations  will  be  accomplished  using  MSOFE  ((4),  (24))  and  MMSOFE. 
As  described  above,  MMSOFE  is  an  adaptation  of  MSOFE  to  enable  MMAE  analysis. 

2.  The  MMSOFE  user  will  already  be  familiar  or  become  familiar  with  the  operation 
of  MSOFE. 

3.  MSOFE  structure  will  be  modified  as  little  as  possible  in  developing  MMSOFE. 

4.  MMSOFE  input  and  output  file  structure,  including  MSOFE  JN,  CTOM,  and  DTOM 
will  be  essentially  identical  to  MSOFE’s  for  ease  of  use. 

5.  MMSOFE  wiU  use  one  MSOFE  JN  input  file  and  produce  one  CTOM  and  one  DTOM 
output  file  for  each  elemental  filter  and  produce  one  CTOM  file  for  the  blended 
MMAE  output  and  parameter  estimation  data. 

6.  MSOFE  simulations  for  elemental  filter  tuning  are  performed  using  covariance  anal¬ 
ysis  and  15-run  Monte  Carlo  analyses  with  statistical  values  averaged  over  the  15 
runs.  Only  Monte  Carlo,  not  covariance,  analysis  wiU  be  used  for  analysis  of  the 
MMAE  results. 

7.  MMSOFE  will  be  structured  to  allow  implementation  of  MMAE,  MMAE-based  con¬ 
trol,  MMAC,  or  DKF  with  minimal  modifications  for  follow-on  research. 

8.  The  MatriXx  (10)  software  will  be  used  only  for  plotting  and  not  for  data  post¬ 
processing;  thus  the  MMSOFE  tool  wiU  be  generally  applicable,  allowing  for  full 
processing  and  plotting  in  computer  environments  that  do  not  include  MatriXx  soft¬ 
ware. 

1.5  Literature  Review 

This  section  contains  a  brief  review  of  literature  pertaining  to  several  FDI  techniques 
other  than  MMAE.  A  review  of  literature  pertaining  to  KF  and  MMAE  theory  and  appli¬ 
cations  is  also  included,  with  references  mainly  from  Maybeck  (18, 19,  20).  The  literature 
review  conclusion  will  discuss  the  usefulness  of  MMSOFE  with  respect  to  the  MMAE-based 
error  detection. 


1-6 


1.5.1  Other  Error  Detection  Methods.  Oae  of  the  simplest  and  most  reliable  failure 
detection  techniques  is  the  use  of  redundant  sensing  elements  for  voting.  Given  a  system 
with  triple  redundancy,  an  algorithm  can  be  easily  written  that  will  compare  the  outputs 
of  each  element,  allowing  them  to  vote  on  the  condition  of  the  other  elements.  Simply 
stated,  if  two  elements  agree  but  a  third  element  totally  disagrees,  the  third  element  is 
considered  faulty  and  is  removed  from  the  system.  Without  the  third  element,  failures 
can  no  longer  be  isolated  by  this  algorithm.  If  the  two  remaining  elements  disagree,  a 
failure  has  been  detected  but  not  isolated,  which  requires  that  both  elements  be  removed 
from  the  system  to  avoid  degradation  of  the  system  performance.  Major  disadvantages  of 
direct  redundancy  lie  in  the  expense,  system  complexity  and  logistics  of  systems  with  the 
requisite  redundant  hardware  (6:pp.  1-2). 

A  more  sophisticated  approach  uses  analytic  redundancy.  By  using  the  physical  and 
dynamic  relationships  between  instruments  on  an  aircraft,  it  is  possible  to  generate  mul¬ 
tiple  sources  of  the  same  information  mathematically.  These  sources  are  used  by  the  FDI 
structure  through  voting  techniques.  Analytic  redundancy  avoids  the  need  for  excessive 
identical  hardware,  by  allowing  non-identical  instruments  to  provide  the  sources  of  sig¬ 
nals  for  the  voting  scheme.  Therefore,  only  the  instruments  normally  required  for  specific 
missions  are  used  without  adding  duplicate  copies  of  instruments.  The  surplus  of  infor¬ 
mation  is  also  used  to  provide  isolation  of  the  failure  without  the  need  for  triple  han-dware 
redundancy  (6:pp.  3-7). 

In  situations  in  which  direct  redundancy  is  impractical,  multiple  Generalized  Like¬ 
lihood  Ratio  (GLR)  testing  (36,  38)  is  an  alternative  which  uses  a  KF  to  compensate  for 
system  failures,  resulting  in  accurate  state  estimates.  Figure  1.2  shows  how  the  sensor, 
the  KF,  and  the  FDI  block  interact.  Three  hypotheses  are  depicted,  with  Hq,  Hi,  and 
representing  no  failure,  jump  failure  and  ramp  failure,  respectively.  The  KF  is  designed 
based  on  the  Hq  hypothesis,  and  the  additional  two  matching  filters  are  designed  based  on 
the  Hi  and  H2  hypotheses  (35).  When  considering  hypotheses  Hi  and  H2,  the  matching 
filters’  design  parameters  determine  what  type  of  failure  is  being  matched  without  prior 
knowledge  of  the  jump  magnitude  or  ramp  slope,  either  or  both  of  which  are  estimated  by 
the  GLR  algorithm.  Each  of  the  matching  filters  monitor  the  KF  residuals  and  computes 
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Figure  1.2.  Multiple  GLR  Testing 

a  maumum  likelihood  estimate  (MLE)  of  the  correctness  of  the  error  hypothesis  mod¬ 
elled  by  that  matching  filter.  The  maximum  likelihood  estimates  (MLE)  determined  by 
each  of  the  hypothesized  matching  filters  are  then  used  in  comparison  test  logic  to  decide 
which  matching  filter  is  correct.  Feedback  to  the  KF  can  also  enable  failure  adaptation. 
Frequently,  a  compensating  bias  as  the  feedback  signal  is  all  that  is  needed  to  continue 
using  the  faulty  sensor.  One  advants^e  of  the  GLR  test  over  some  other  FDI  algorithms  is 
that  prior  total  knowledge  of  the  failure  magnitude  is  unnecessary.  GLRs  and  MLE’s  are 
presented  in  detail  by  Willsky  and  Jones  (36,  37,  38). 

The  simple  hypotheses  for  fail/no-fail  conditions  may  not  provide  the  robustness 
needed  to  detect  and  adapt  to  certain  failures.  Detection  of  ramp  fmlures  is  particularly 
difficult  when  the  ramp  rate  varies  significantly.  A  time  delay  also  exists  in  detecting  ramp 
failures.  This  delay  occurs  because  the  deceiving  signal  is  slowly  moving  away  from  the 
desired  signal  and  takes  longer  to  cross  the  failure  threshold  than  with  a  jump  failure. 
These  failure  threshold  checks  might  be  as  simple  as  monitoring  the  residuals  to  see  when 
they  exceed  some  threshold  value  or  they  could  be  more  complex  methods.  The  residuals 
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are  particularly  succeptible  to  change  when  failures  are  induced  as  biases  or  ramps  directly 
in  the  measurement  or  in  the  measurement  noise  variance. 

Considering  a  filter  which  hypothesizes  no  failure,  in  the  case  of  a  true  measurement 
bias/ramp,  the  residual  itself  will  be  consistently  offset  from  its  prior  value  (circa  zero 
hopefully)  and  in  the  case  of  the  measurement  noise  variance  bias/ramp  failure,  the  actual 
covariance  of  the  observed  residuals  increases  when  the  failure  is  present.  By  adding  filters 
designed  to  match  the  ramp  failure,  these  longer  delays  can  possibly  be  avoided.  It  is 
also  desirable  to  design  the  filter  based  on  a  specific  assumed  time  of  the  failure.  If  the 
filter  were  customized  to  look  for  a  failure  at  a  specific  instant  in  time,  then  it  would 
have  a  better  chance  of  accurate  detection,  but  this  idea  results  in  a  very  large  bank  of 
filters  for  long  periods  of  time.  Additionally,  GLR  could  estimate  the  time  of  onset  of  a 
ramp  as  well  as  estimate  the  slope  of  the  ramp.  The  disadvantage  of  adding  filters  is  the 
increase  in  computations  required  for  multiple  filters.  A  solution  to  this  problem  uses  a 
set  number  of  filters  over  a  window  of  time  in  which  each  filter  differs  slightly,  so  that  one 
filter  will  hopefully  match  the  real  world.  The  number  of  filters  remains  constant,  and  the 
performance  of  the  FDI  system  is  often  maintained  (38:pp.  606-608).  The  justification  for 
using  sliding  windows  is  discussed  by  Willsky  (38:pp.  604-605). 

Another  FDI  method  based  on  residual  monitoring  is  a  chi-square  test,  which  is 
similar  to  GLR  testing  in  that  it  calculates  a  xi^t)  random  variable  based  on  the  filter 
residuals,  given  by 

x(l.)=  i:  p’'(l,)A-'((,)r(<,)  (1.1) 

with  N  being  the  size  of  a  sliding  window,  r  being  the  vector  of  the  residuals,  and  A{tj) 
being  the  covariance  of  the  residuals.  One  difference  between  these  two  algorithms  is  that 
GLR  tests  are  explicitly  functions  of  the  assumed  system  dynamics  model  and  chi-square 
tests  are  also  but  only  indirectly.  In  any  case,  the  chi-square  test  can  only  determine  if  a 
given  hypothesis  is  correct.  A  second  major  difference  is  that  chi-square  tests  do  not  try 
to  match  the  failure  and  only  have  one  hypothesis,  making  it  a  binary  test  for  fail/no-fail. 
Without  the  ability  to  distinguish  between  types  of  failures,  the  chi-square  test  is  not 
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Figure  1.3.  One  Possible  Configuration  of  Multiple  Model  Adaptive  Estimation 


useful  for  isolation.  However,  cbi-square  testing  is  easy  to  implement,  runs  quickly,  and 
can  provide  the  first  level  of  failure  detection  in  a  multi-level  FDI  scheme. 

1.5.2  Multiple  Model  Adaptive  Estimation  (MMAE).  Of  primary  interest  for  this 
research  is  this  final  technique,  the  use  of  multiple  models  that  represent  the  dynamics  of 
the  system  to  include  various  failure  conditions.  This  technique  is  somewhat  analogous  to 
multiple  GLR  testing  in  that  it  explicitly  accounts  for  multiple  hypotheses,  but  it  differs  in 
its  structure  and  decision  making  process  and  in  its  capability  of  improved  performance. 
Unlike  the  multiple  GLR  testing  which  modeled  a  variety  of  failures  using  matching  filters, 
the  MMAE  models  the  dynamic  nature  of  the  system  and  its  sensors  to  represent  their 
behavior  in  the  presence  of  a  set  of  particular  assumed  failures.  Figure  1.3  shows  elemental 
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Kalman  filters  which  are  designed  to  model  the  various  failure  conditions.  The  elemental 
filter  residuals  are  used  to  determine  which  filter  best  models  the  aircraft  and  its  sensors 
at  the  current  time  and  thus  to  indicate  which  hypothized  failure  has  occurred.  A  jump 
or  ramp  might  actually  be  modeled  by  several  elemental  filters,  each  of  which  models  a 
different  magnitude  of  bias.  It  should  be  noted  that  ramp  failures  probably  don’t  require 
rsimp  hypotheses,  but  given  bias  hypotheses,  the  MMAE  algorithm  can  switch  sequentially 
from  one  assumed  bias  value  to  another.  A  probability  of  adequacy  of  the  model  assumed 
in  each  elemental  filter  ranging  from  zero  to  one  is  computed  for  each  elemental  filter 
and  is  used  to  weight  that  filter’s  state  estimates  to  form  a  probability-weighted  average 
as  an  output  estimate.  A  Bayesian  form  of  MMAE  filter  blends  (weighted-averages)  the 
probability-weighted  estimates  of  the  respective  elemental  filters  to  produce  the  optimal 
blended  estimates.  A  Maximum  A-Posteriori  (MAP)  form  of  MMAE  algorithm,  in  con¬ 
trast,  would  select  the  one  elemental  filter  with  the  highest  probability,  and  output  that 
single  filter’s  estimate  instead.  Similarly  the  parameters  which  vary  between  elemental 
filters  can  be  weighted  with  the  corresponding  probabilities  to  estimate  the  true  value  of 
that  parameter.  This  MMAE  blending  approach  allows  for  partial  sensor  failures  or  com¬ 
binations  of  failure  types,  as  in  the  blending  together  a  non-failure  filter  and  a  filter  which 
assumes  a  jump  bias  of  magnitude  J  in  a  specific  sensor.  This  blending  could  handle  an 
ctctual  jump  value  between  zero  and  J.  However,  the  discretization  of  values  of  J  cannot  be 
too  coarse  or  else  no  filter  might  be  close  enough  to  the  real-world  bias  to  generate  well- 
behaved  residuals.  Probabilities  of  accuracy  ranging  from  zero  to  one  correspond  to  filter 
models  that  are  deemed  to  be  completely  inaccurate  or  perfectly  accurate,  respectively. 
A  state  estimate’s  probability  corresponds  directly  to  the  decimal  portion  of  that  state 
estimate  which  contributes  to  the  blended  estimate.  MMAE  is  described  in  detail  by  May- 
beck  (20)  and  has  been  used  for  numerous  applications  as  discussed  in  Section  1.5.3.  An 
important  adaptation  of  MMAE  is  that  of  either  MMAE-based  control  or  multiple  model 
adaptive  control  (MMAC)  for  system  stability  and  failure  correction  (related  diagrams 
are  shown  in  Chapter  II).  MMSOFE  should  be  easy  to  modify  with  a  user  subroutine  to 
accommodate  MMAE-based  control  or  MMAC. 
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1.5.3  MMAE  Applications.  As  mentioned  previously  in  Section  1.1,  MMAE  has 
been  used  for  many  different  applications  and  most  notably  at  the  Air  Force  Institute  of 
Technology  (AFIT)  (7,  8,  11,  21,  26,  29,  30,  33),  Alphatech  Inc.  (13),  and  The  Computer 
Technology  Institute  (Patras,  Greece)  (14).  MMAE,  MMAE-based  control  and/or  MMAC 
have  been  the  method  of  choice  for  a  number  of  applications-oriented-research  projects: 

1.  Detect  and  identify  sensor  and  aircraft  flight  control  sensor  and/or  control  surface 
failures  (8,  21). 

2.  Reconfigure  an  aircraft  flight  control  system  when  system  failures  occur  (33). 

3.  Select  among  different  first-  and  second-order  target  acceleration  models  for  use  in  a 
correlator/forward-looking-infrared  tracker  for  airborne  targets  (26). 

4.  Implement  an  adaptive  model-based  control  for  a  robotic  arm  (29). 

5.  Study  target  tracking  and  weapon  lead  prediction  using  MMAE  (13). 

6.  Implement  a  moving-bank  subset  of  multiple  models  for  flexible  spacestructure  con¬ 
trol  in  the  face  of  uncertain  parameters  to  describe  the  bending  modes  (7,  11,  30). 

1.5.4  Literature  Review  Conclusion.  The  simple  error  detection  methods  presented 
in  Section  1.5.1  might  be  effective  for  detecting  some  failures,  but  may  be  inadequate  for 
more  complex  systems  such  as  detecting  errors  in  GPS  pseudorange  in  a  GPS-aided  inertial 
navigation  system,  due  to  jamming  or  spoofing.  These  errors  are  small  in  comparison  to 
some  of  the  other  GPS  errors  and  are  therefore  difficult  to  detect.  While  such  errors  would 
degrade  the  GPS  performance,  the  cause  of  the  failures  would  still  be  undetectable  with 
such  a  simple  error  detection  method  (5). 

Direct  redundancy  techniques  using  multiple  GPS  receivers  or  extra  satellites  to 
provide  redundancy  are  expensive  and  normally  impractical,  but  not  impossible.  However, 
analytical  redundancy  is  inherent  to  integrated  navigation  systems  in  that  the  sensors  (i.e., 
GPS,  INS,  barometric  altimeter,  and  etc.)  can  be  used  to  provide  multiple  sources  of  the 
same  information. 

Multiple  GLR  testing  and  MMAE  are  prime  candidates  for  an  EDI  system.  Both  of 
these  techniques  possess  the  EDI  versatility  needed  for  complex  systems,  but  their  major 
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difference  is  in  the  type  of  multiple  filters  required  by  each  and  the  accuracy  of  the  state 
estimates.  Multiple  GLR  testing  uses  a  bank  of  matching  filters  that  are  somewhat  less 
complicated  than  the  bank  of  KFs  used  in  MMAE.  Quicker  response  might  be  anticipated  in 
MMAE  because  the  elemental  filters  are  running  in  parallel  at  all  times  based  on  competing 
hypotheses  versus  only  one  filter.  Vasquez  (36)  considered  both  techniques  as  candidates  for 
his  EDI  algorithm  and  chose  to  use  the  GLR  tests  and  chi-square  tests.  He  used  the  GLR 
algorithm  only  for  failure  isolation,  but  MMAE  should  be  capable  of  simultaneous  failure 
detection,  failure  isolation,  and  generation  of  blended  estimates  which  are  superior  to  the 
estimates  from  any  of  the  MMAE  elemental  filters  or  those  resulting  from  the  GLR  method. 
Performance  evaluation  of  those  other  approaches  can  be,  and  has  been,  established  with 
the  standard  MSOFE  (or  similar)  analysis  tool.  Therefore,  the  MMSOFE-MMAE  tool 
developed  in  this  thesis  will  make  it  much  easier  to  compare  MMAE  FDI  results  with 
results  from  any  of  the  aforementioned  methods. 

1.6  Methodology 

The  main  steps  of  this  research  are  explained  in  Sections  1.6.1  through  1.6.5.  Steps 
one  and  two  will  satisfy  the  primary  goal  of  the  thesis,  while  step  three  is  aimed  at  the 
follow-on  goal  included  as  a  recommendation  for  future  research. 

1.6.1  MMAE  Implementation  Options.  There  were  several  methods  considered  for 
implementing  MMAE  using  MSOFE. 

1.  Generating  new  general-use  MMAE-specific  software  would  always  be  an  option, 
either  with  or  without  using  segments  of  MSOFE.  This  would  mean  the  generation 
of  new  software  which  would  not  be  familiar  to  MSOFE  users  and  therefore  require  a 
greater  learning  curve  for  a  new  user  than  would  be  required  with  an  MSOFE-similar 
approach  like  in  4  below. 

2.  One  obvious  method  is  to  run  MSOFE  serially,  once  for  each  elemental  filter  and  then 
use  post-processing  to  calculate  the  probabilities,  blended  state  estimates,  parame¬ 
ter  estimates  and  other  desired  values.  For  simple  MMAE,  this  is  a  viable  option, 
although  the  post-processing  algorithms  would  probably  be  re-invehted  for  each  ap- 
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plication  because  the  particular  software  package  used  for  post-processing  might  not 
be  available  or  familiar  to  the  user.  Morever,  true  feedback  for  MMAE-based  control 
would  still  not  be  possible. 

3.  Another  method  using  MSOFE  would  be  to  add  a  set  of  states  for  each  elemental 
filter  to  the  “filter”  model  used  in  MSOFE,  thus  forming  one  large  filter.  This 
one  filter  would  be  equivalent  to  the  bank  of  elemental  filters,  and  the  blending 
computation  would  have  to  be  augmented  to  that  “single  filter”  algorithm.  With 
this  implementation,  the  number  of  filter  states  in  the  MSOFE  model  would  equal 
the  number  of  states  required  by  the  problem  multiplied  by  the  number  of  elemental 
filters.  This  approach  would  be  possible  but  the  total  number  of  states  would  grow 
rapidly  for  more  complex  problems.  Further,  the  MMAE  probability  computing  and 
estimate  blending  algorithms  would  need  to  be  programmed  separately  or  possibly 
as  even  more  additional  states.  Therefore,  a  four-state  MMAE  problem  with  four 
elemental  filters  might  require  twenty  states.  MMAE-based  control  could  be  possible 
but  implementation  would  be  complex. 

4.  Last  is  a  method  which  would  truly  implement  MMAE.  It  would  calculate  a  prop¬ 
agation  and  update  cycle  for  each  elemental  filter,  followed  by  the  MMAE-specific 
algorithms,  before  any  filter  would  begin  another  propagation  cycle.  This  method 
was  chosen,  because  MSOFE  could  be  modified  to  preserve  its  features  and  because 
this  method  allows  any  multiple  model  application,  be  it  MMAE,  MMAE-based 
control,  or  DKF  to  be  run  with  minimal  user  modifications. 

1.6.2  Software  Development.  As  mentioned  in  the  literature  review,  multiple  GLR 
testing  has  an  advantage  over  MMAE  of  needing  only  one  KF,  while  MMAE  requires  a 
bank  of  such  filters  running  in  parallel.  Contrarily,  MMAE  has  the  advantages  of  parameter 
estimation,  FDI,  improved  state  estimates  and  control  in  the  MMAE-based  control  form. 
Morever,  the  very  existence  of  a  parallel  bank  of  filters,  each  based  on  a  different  hypothesis 
and  each  producing  residuals  associated  with  a  particular  hypothesis,  inherently  aids  the 
hypothesis  testing  process  in  the  MMAE.  Many  previous  AFIT  Kalman  filter  theses  utilized 
single-filter  MSOFE  runs  on  the  VAX,  but  with  MMSOFE  running  multiple  filters  in 
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parallel,  MMAE  on  the  VAX  would  take  entirely  too  long.  Therefore,  the  satellite  orbit 
problem  was  used  with  MSOFE  and  MMSOFE  on  the  SPARC  workstations.  MSOFE  was 
modified  into  MMSOFE  to  run  MMAE,  and  MMSOFE  was  verified  using  the  satellite 
orbit  problem  as  presented  by  Maybeck  (19:pp.  46-48)  and  in  the  MSOFE  users’  manual 
(24)  as  a  test  case.  Further,  the  data  output  from  the  each  elemental  filter  was  compared 
with  the  output  for  the  corresponding  filter  run  with  MSOFE.  This  permitted  a  bit-by-bit 
comparison  of  the  data  produced  by  MMSOFE  with  that  of  MSOFE,  demonstrating  that 
the  implementation  and  performance  assessment  of  the  elemental  filters  within  the  MMAE 
were  error-free. 

1.6.3  Preliminary  Studies.  The  satellite  orbit  problem  used  as  the  test  case  for  this 
research  has  range  and  time  measured  in  generic  radius-units  and  time-units.  This  study 
covers  both  ramp  and  jump  failures  in  the  range  measurement  and  range  measurement 
noise  variance  at  different  magnitudes  and  rates,  as  shown  in  Table  1.1.  The  primary 
goal  of  this  satellite  orbit  problem  will  be  software  validation  with  the  intrinsic  goals  of 
accomplishing  FDI  in  the  range  estimate  of  the  satellite.  No  effort  wiU  be  made  to  adapt 
to  any  failures. 


Table  1.1.  Proposed  Orbit  Problem  Failure  Cases 


Error 
Case  # 

Failure  Type 

In  Truth  Model 

Failure 

Induced 

Elemental  Filter 

1 

Step  in 

Measurement  Noise 

Rs  =  Rso  +  ^Rso 

T  >  25 

i2/o(l  +  RIA-Siij) 

2 

Ramp  in 

Measurement  Noise 

25  <  T  <  35 

r  >  35 

R/o(l  +  BIASrj) 

3 

Step/Return  in 
Measurement  Noise 

Rs  —  Rso 

Rs  =  Rso  +  ^Rso 

T  <  25,  r  >  35 
25  <  r  <  35 

R/o(l  +  B/A5fi,) 

4 

Step  in 
Measurement 

r  >  10 

Z,  -t-  BIASz, 

5 

Ramp  in 
Measurement 

10  <  T  <  20 

T  >  20 

Zf  -1-  BlASz, 

To  explain  the  entries  in  Table  1.1  further: 


T  =  Current  simulation  time. 
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Rs  =  Truth  model  (real-world  or  “System”)  range  measurement  noise  variance. 

RgQ  =  Nominal  value  of  Rg. 

Rfo  =  Elemental  filter  nominal  range  measurement  noise  variance. 

Zs  =  Truth  model  range  measurement. 

Z/  =  Elemental  filter-predicted  range  measurement. 

BIASrj  =  Elemental  filter  scale-factor  on  Rfo  for  bias  added  to  Rjo  (BIASrj  = 
{6, 3,0, -.5, -.75}). 

BIASz,  =  Elemental  filter  biases  added  to  Z/o 
{BIASz,  =  {.1,. 05,0, -.05,-.!}). 

The  onset  and  relaxation  times  for  the  failures  were  chosen  so  there  would  be  enough 
time  both  before  and  after  the  failures  for  the  elemental  filters’  residuals  to  stabilize  at 
some  “steady  state”  value,  thus  resulting  in  MMAE-calculated  probabilities  shifting  from 
baseline  to  the  appropriate  elemental  filters  after  induction  of  the  failure.  The  magnitude  of 
the  failures  in  the  truth  models  were  chosen  to  be  large  enough  to  make  it  possible  for  good 
discrete  failure  detection  and  isolation  among  the  elemental  filters,  while  still  not  being 
so  large  as  to  be  totally  unrealistic.  The  ms^nitudes  chosen  for  BIASr^  and  BIASz, 
were  chosen  to  make  the  failure  detection  obvious,  with  the  largest  values  in  each  case 
matching  the  magnitude  of  the  truth  model  failure  bias  or,  in  the  case  of  ramps,  matching 
the  final  (sustained)  magnitude.  The  other  magnitudes  of  BIASr^  were  chosen  such  that 
the  resulting  noise  variance  would  be  |  of  the  prior  value.  Once  again  for  BIASzf,  the 
largest  magnitude  matches  that  in  the  system  model,  with  the  second  value  being  |  the 
first  value,  followed  by  the  baseline  filter  which  hypothesizes  no  failure.  The  last  two  are 
just  the  negative  (mirror  image)  of  the  first  two. 

1.6.4  Software  Validation  Tests.  Because  development  of  the  MMAE  software  tool, 
MMSOFE,  is  the  purpose  of  this  thesis,  software  validation  is  paramount  to  this  research, 
particularly  since  MMSOFE  will  be  an  entirely  new  adaptation  of  MSOFE  which  is  not 
without  risk.  The  steps  in  testing  the  MMSOFE  software  include: 
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1.  Modify  the  orbit  MSOFE  SPARC  code  carrently  available  for  the  chosen  simulated 
errors  (see  Table  1.1). 

2.  Conduct  the  MSOFE  runs  shown  in  Table  1.1  and  evaluate  the  results. 

3.  Tune  the  filters  which  represent  each  failure  type  via  MSOFE  covariance  analysis 
runs. 

4.  Generate  the  new  SPARC  MMSOFE  orbit  code  by  modifying  the  SPARC  MSOFE 
orbit  code  for  the  MMAE  algorithm. 

5.  Determine  the  filter  models  and  elemental  filters  required  in  MMAE  for  the  failure 
types  shown  in  Table  1.1. 

6.  Partially  validate  the  MMSOFE  software  by  using  multiple,  but  identical  elemental 
filters  with  the  MMAE  probability  and  blending  computations  disabled.  Compare 
the  output  with  that  obtained  from  MSOFE. 

7.  Validate  the  MMSOFE  softwaire  by  using  multiple,  but  distinct  elemental  filters  and 
compare  the  elemental  filter  results  to  data  obtained  from  matching  MSOFE  runs. 

8.  Generate  MMSOFE  to  employ  Bayesian-blending  for  the  optimal  state  and  param¬ 
eter  estimation. 

9.  Run  MMSOFE  with  no  Hohmann  orbit  transfer  for  the  first  four  failure  types  shown 
in  Table  1.1  and  compare  the  blended  state  estimates  to  the  output  of  the  elemental 
filters. 

10.  Possibly  run  MMSOFE  with  Hohmann  orbit  transfer  for  the  first  four  failure  types 
shown  in  Table  1.1  and  compare  the  blended  state  estimates  to  the  output  of  the 
elemental  filters. 

1.6.5  Follow-on  Goals.  As  (L  cussed  above,  one  follow-on  goal  is  to  adapt  MM- 
SOFE  to  provide  the  optimal  corrective  feedback  and  thus  make  the  MMAE  system  into 
a  truly  adaptive  MMAE-based  control  system  or  MMAC.  Another  follow-on  goal  is  to 
investigate  adapting  MMSOFE  for  Distributed  Kalman  Filter  (DKF)  applications  as  well. 
Further,  MMSOFE  could  be  used  merely  to  provide  multiple  tuning  runs  by  using  elemen¬ 
tal  filter  with  identical  models  but  different  tuning  values  versus  making  multiple  MSOFE 
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tuning  runs  with  the  different  tuning  values.  An  added  benefit  to  using  MMSOFE  in  this 
manner  would  be  the  parameter  estimation  to  aid  in  determining  the  correct  tuning  values. 

Lastly,  a  foUow-on  application-oriented  goal  is  to  use  MMSOFE  on  the  Navigation 
Reference  System  (NRS)  FDI  problem  investigated  by  Vasquez  (36)  who  used  GLR  and 
chi-square  techniques.  Adaptive  corrective  feedback  could  improve  the  state  estimates  and 
FDI  capabilities  of  a  MMAE-blended  NRS  system. 

1.6.6  Stopping  Criteria.  The  primary  stopping  criteria  correlate  to  the  steps  in 
Section  1.6.4.  The  stopping  criteria  for  software  development  and  validation  are  as  follows; 

1.  MSOFE  orbit  estimation  problem  tuning  runs  with  the  errors  shown  in  Table  1.1 
will  be  complete  when  the  true  error  standard  deviations  computed  in  covariance 
analysis  roughly  match  the  filter-predicted  values. 

2.  MMSOFE  algorithm  wiU  be  validated  when  the  output  for  each  elemental  filter 
exactly  matches  that  produced  by  running  the  identical  model  with  MSOFE,  first 
when  all  elemental  filters  match  and  then  when  they  differ. 

3.  The  MMAE  blending  algorithm  will  be  complete  when  the  blended  state  estimates 
are  reasonable  in  comparison  to  the  state  estimates  produced  by  the  elemental  filters. 

1.7  Overview  of  Thesis 

Chapter  II  presents  the  detailed  theory  used  in  the  research.  KF  theory  is  discussed, 
including  time-discretization  of  the  dynamics  for  sampled  data  KF  implementation  and 
with  special  attention  on  MMAE,  MMAE-based  control  and  MMAC,  and  some  attention  to 
DKF.  The  basics  of  the  FDI  algorithms  are  discussed,  including  the  equations  implemented 
for  the  MMAE. 

Chapter  III  contains  the  detailed  description  of  the  MMSOFE  adaptation  of  MSOFE, 
to  include  the  theoretical  basis,  the  software  development  and  detailed  modifications  to 
MSOFE.  The  basics  of  the  orbit  estimation  problem  are  discussed,  as  well  as  the  results 
from  the  MMSOFE  preliminary  verification  simulations. 
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Chapter  IV  describes  results  of  the  MMSOFE  orbit  problem  simulations.  Analysis 
of  the  elemental  filter  performance,  the  blended  state  estimates,  the  probabilities  and 
parameter  estimation,  and  the  FDI  results  will  also  be  discussed. 

Chapter  V  summarizes  the  research  and  presents  conclusions  and  recommendations 
with  regard  to  using  MMSOFE  for  other  multiple  model  applications. 
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II.  Kalman  Filtering,  MMAE,  and  Failure  Detection 


2.1  Overview 

This  chapter  presents  the  fundamental  theory  of  Kalman  filtering,  both  when  used 
as  a  single  filter  and  as  employed  in  a  bank  of  elemental  filters  in  the  MMAE  context. 
Both  MSOFE  and  MMSOFE  can  run  linear  as  well  as  extended  Kalman  filter  simulation. 
Because  the  satellite  orbit  problem  used  for  MMSOFE  development  in  this  thesis  uses  an 
extended  Kalman  filter  (EKF)  algorithm,  a  suitable  discretization  method  of  the  system 
dynamics  must  be  used.  The  method  of  discretizing  the  filter  equations  for  continuous-time 
dynamics  will  be  described.  MMAE-related  FDI  methods  and  MMAE  blended  estimation 
will  be  discussed.  Stochastic  adaptive  contollers  in  the  form  of  MMAE-based  controllers 
and  Multiple  Model  Adaptive  Controllers,  and  Distributed  Kalman  Filtering  will  also  be 
discussed  with  respect  to  MMSOFE.  Finally,  filter  tuning  and  selection  of  lower  bounds 
for  MMAE-computed  hypothesis  conditional  probabilities  will  be  discussed  briefly. 

2.2  The  Extended  Kalman  Filter 

Because  of  the  dynamic  nature  of  the  chosen  problem  of  interest  and  for  the  sake  of 
more  general  applicability,  an  Extended  Kalman  Filter  (EKF)  was  chosen  to  provide  state 
estimates.  The  EKF  allows  for  nonlinear,  time-varying  dynamics  and/or  measurement 
equations.  In  simple  linearized  Kalman  filtering  (LKF),  these  equations  are  linearized 
through  approximation  techniques  about  a  fixed  nominal  trajectory  to  form  the  LKF.  The 
LKF  is  the  conceptual  basis  for  the  EKF,  but  the  EKF  is  continually  relinearized  about 
the  most  recent  state  estimate  rather  than  about  a  fixed  nominal  trajectory. 

2.2.1  The  Sampled  Data  Kalman  Filter.  Let  the  system  model  be  expressed  as  a 
state  equation  of  the  form: 


x(t)  =  f[x(«),t]  +  G(0w(f)  (2.1) 

where  the  state  dynamics  vector  f[x(t),t]  is  a  nonlinear  function  of  the  state  vector  x(f) 
and  time.  It  is  noteworthy  to  point  out  that  this  discussion  of  LKF  and  EKF  equations 
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correlates  very  closely  with  the  discussion  in  Section  9.5  of  Maybeck  (19).  However,  there 
is  one  simplification  in  the  state  dynamics  vector,  f[x(t),u(t),t]  Maybeck’s  Equation  (9- 
34b),  where  u(t)  represents  the  deterministic  control  inputs.  In  general  extended  Kalman 
filtering  and  with  MSOFE  and  MMSOFE,  the  u(t)  inputs  are  non-zero  and  non-constant. 
For  the  MMSOFE  development  simulations,  however,  it  is  assumed  that  u(t)  =  0,  which 
corresponds  to  no  deterministic  control  inputs.  With  that  aside,  let  the  process  noise,  w(t) 
be  a  white  Gaussian  noise  with  mean: 


E  {w(t)}  =  0 


(2.2) 


and  noise  strength  Q(t)  defined  by: 

E  {w(t)w^(t  +  r)}  =  Q(t)«(T)  (2.3) 

The  discrete-time  measurements,  z(t,  )  are  modeled  as  a  nonlinear  vector  of  functions 
of  the  state  vector  and  time,  h[x(t,),t,]  plus  additive  white  Gaussian  noise; 


z(<<)  =  li[x(t<),<<]  -f  v(t<)  (2.4) 

where  v(ti)  is  a  zero-mean  discrete-time  white  Gaussian  noise  of  covariance  R(t<)  defined 
by: 

E  {v(t<)v^(t^)}  = 

and  h[x(ti),tj]  is  the  nonlinear  observation  (measurement)  vector.  The  LKF  is  based 
on  perturbation  states  about  a  fixed  nominal  state  trajectory  x„(t)  satisfying  the  initial 
condition,  x„(to)  =  Xno  (some  value,  assumed  to  be  known),  and 

x„(0  =  f(x„(0,^]  (2.6) 

where  f[-,  •]  is  shown  in  Equation  (2.1 ).  The  corresponding  nominal  discrete-time  measure¬ 
ments  for  the  nominal  trajectory  are  based  on  the  nominal  states  and  time,  but  without 


R(t,)  for  ti  =  tj 
0  for  ti  /  tj 


(2.5) 
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aay  additive  noise: 


*»(<<)  =  h[x„(t,),tj] 


(2.7) 


The  perturbation  states  are  found  by  subtracting  the  nominal  states  in  Equation  (2.6) 
from  the  original  states  in  Equation  (2.1): 

[x(0  -  Xnit)]  =  flx(t),  t]  -  flx„(t),  t]  +  G(0w(0  (2.8) 

Equation  (2.8)  is  approximated  to  first  order  through  a  truncated  Taylor  series  expansion: 

Sx(t)  =  F  [t;  x„{t)]  6x(t)  +  G(t)w(f)  (2.9) 

where  ^x(t)  are  the  perturbation  states.  The  definitions  for  G(t)  and  w(f)  are  unchanged 
and  the  new  linearized  dynamics  matrix  F[t;  x„(t)]  is  found  by  taking  partial  derivatives 
of  l[x(f ),  t]  with  respect  to  x{t)  and  evaluated  along  the  nominal  values  for  the  trajectory 
*«(f)5 

F(I1X,(I)|  = 

The  discrete-time  perturbation  measurements  are  similarly  approximated  to  first  order 
from  the  measurement  difference  equation: 

[z(«)  -  *„(f)]  =  h(x(f.),«i]  -  h[x„(t,),fi]  +  y{ti)  (2.11) 

yielding  the  perturbed  form: 

Sz{ti)  =  H[<<;x(f,)]«x„(t)  +  v(t,)  (2.12) 


(2.10) 

X=X,(t) 


The  same  partial  derivative  methods  used  to  derive  the  linearized  state  dynamics  matrix 
F[f ;  x„(f)]  are  used  again  to  derive  the  linearized  observation  matrix: 


H(fj;x„(<<)] 


ah[x(t<),f<] 

dx 


lx=x.(<.) 


(2.13) 
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Because  the  LKF  generates  error  state  estimates,  Sx{t),  they  must  be  added  to  the  nominal 
states  to  provide  whole  states  estimates  x(t)  in  the  form: 

x(0  =  x„(t)  +  Sx{t)  (2.14) 


However,  if  the  nominal  and  the  “true”  trajectories  differ  too  greatly,  a  linearized  Kalman 
filter  can’t  really  be  used  because  the  condition  for  negli^bility  of  higher  order  terms  in  the 
Taylor  series  expansions  has  been  violated.  In  such  situations,  an  extended  Kalman  filter 
can  be  implemented  by  linearizing  about  the  current  state  estimate  x  rather  than  about 
the  nominal  trajectory  x„.  The  following  sampled  data  extended  Kalman  filter  equations 
use  the  notation  t/ti  to  represent  the  value  of  a  ^ven  variable  at  time  t,  conditioned  on  the 
measurements  taken  through  the  time  t,.  Also,  represents  the  value  after  propagation 
but  prior  to  the  measurement  update  at  sample  time  f,-,  and  tf  corresponds  to  the  value 
after  the  measurement  update.  The  state  estimates  x  and  covariance  values  P(t/ti)  are 
propagated  from  U  to  ti+i  by  solving  the  following  differential  equations: 


P(t/<.)  =  F[t;x(«/ti)]P(t/t<)  +  +  G(f)Q(0G’'(0 


where 


F[t;x(t//,)]  = 


dx 


and  initial  conditions  are  given  by: 


(2.15) 

(2.16) 

(2.17) 


=  *(<?■) 

(2.18) 

P((,/I,)  =  P(tt) 

(2.19) 

The  discrete-time  measurements  are  processed  in  the  EKF  through  update  equations: 

K(<.)  =  P(tr)H’’[t<;x(tr)]  {H[ti;x(tr)]P(tr)H’’[fi;x(fr)]  +  R(^.)}"‘  (2.20) 
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x(t+  )  =  x(t" )  +  K(t.)  {*<  -  h[xiti  ), f.]}  (2.21) 

P(it)  =  P(fr)  _  KiU)H[ti-Mt7mt-)  (2.22) 


where 


H[<.;x(«r)l 


ah[x(<<),<.] 


dx 


lx=*(»") 


and  K(<j)  is  the  discrete-time  Kalman  filter  gain. 


(2.23) 


It  is  noteworthy  that  the  covariance  updates,  Equation  (2.22),  are  actually  imple¬ 
mented  in  MSOFE  using  U  —  D  covariance  factorization  (24).  This  covariance  update 
form  is  shown  by: 


p(i.-)  =  u(tr)D(i-)tj-'(ir) 

(2.24) 

p(tr) = 

(2.25) 

in  which  the  U  is  a  unitary  (diagonal  terms,  which  are  the  eigenvalues,  equal  one),  upper 
triangular  matrix  and  the  D  matrix  is  dis^onal.  The  numerical  advantages  inherent  to 
U  —  D  covariance  factorization  are:  guaranteed  nonnegativity  of  the  covariance  matrix, 
numerical  accuracy,  and  numerical  stability  (18).  Because  MMSOFE  is  an  MSOFE  deriva¬ 
tive,  the  same  U  —  D  covariance  factorization  advantages  are  preserved  in  MMSOFE. 

2.3  Multiple  Model  Adaptive  Estimation 

MMAE  is  capable  of  failure  detection  and  isolation,  optimal  state  estimation  and 
parameter  estimation.  MMAE  is  composed  of  multiple  Kalman  filters  running  in  paral¬ 
lel  (simultaneously),  as  discussed  in  Section  1.5.2.  Each  elemental  filter  models  one  of  the 
possible  configurations  of  the  systems.  These  configurations  may  be  intentional,  such  as  el¬ 
emental  filters  based  upon  different  possible  system  parameter  values,  or  the  configurations 
may  represent  some  type  of  failure,  such  as  those  caused  by  GPS  jamming  or  spoofing.  If 
these  failures  are  of  unknown  magnitudes,  the  MMAE  must  estimate  both  the  states  and 
the  failure  magnitude  by  using  elemental  filters  with  different  assumed  failure  magnitudes. 
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This  is  only  slightly  different  from  estimation  of  an  unknown  system  parameter.  Given 
jumps  or  ramps  in  the  measurement  signals  or  in  the  measurement  noise  variances  as  fail¬ 
ure  types,  the  elemental  filters  will  use  different  assumed  bias  values,  and  the  magnitude 
of  the  true  failure  can  be  estimated  by  merely  summing  the  probability-weighted  assumed 
values  that  were  used  in  the  elemental  filters.  There  is  no  need  to  hypothesize  elemental 
filters  based  on  a  ramp  failure  of  an  assumed  slope  and  assumed  time  of  onset;  as  the  true 
ramp  occurs,  the  elemental  filter  with  the  bias  value  that  most  closely  matches  the  current 
ramp  value  ought  to  achieve  the  highest  probability  weighting.  This  magnitude  estimation 
develops  naturally  after  calculating  the  conditional  probabilities  for  each  elemental  filter. 
These  conditional  probabilities,  also  sometimes  referred  to  as  hypothesis  conditional  prob¬ 
abilities,  will  be  used  as  weights  for  blending  the  state  estimates  from  the  corresponding 
elemental  filters  to  determine  the  blended  MMAE  estimates.  These  conditional  probabil¬ 
ities  are  also  used  to  estimate  the  uncertsdn  parameter  value,  which  in  the  case  of  FDI  is 
the  true  magnitude  of  the  jump  or  ramp. 

Out  of  all  of  the  possible  filters  existing  in  the  entire  parameter  space,  the  approach 
in  this  research  will  be  to  use  several  elemental  Kalman  filters  which  differ  only  by  the  as¬ 
sumed  measurement  noise  variance  or  measurement  signal  bias  magnitude.  Measurement 
signal  biases  would  model  measurement  jumps  and/or  ramps  as  might  result  from  inten¬ 
tional  spoofing  by  an  adversary,  and  increased  measurement  noise  variance  woiild  model 
jamming-type  real  world  measurement  noise.  Nominally,  for  the  case  of  measurement  sig¬ 
nal  biases,  two  filters  with  different  positive  bias  magnitudes,  a  matched  set  with  negative 
magnitudes,  and  one  filter  with  bias  magnitude  =  0  can  be  assumed  (i.e.,  K  =  5,  where  K 
is  the  number  of  elemental  Kalman  filters).  Using  varying  variances  of  measurement  noise 
instead  of  biases  in  another  set  of  elemental  filters  could  similarly  model  jamming  noise. 
Thus,  the  probabilities  determined  for  each  of  these  filters  will  correspond  to  their  relative 
closeness  to  the  truth  model  (real  world)  characteristics.  The  conditional  probabilities 
corresponding  to  a  good  filter  model  (small  residuals)  will  be  greater  than  the  conditional 
probabilities  corresponding  to  a  poor  filter  model  (large  residuals).  Thus  the  elemental 
filter  state  estimates  are  weighted  in  proportion  to  the  probability  that  the  specific  ele¬ 
mental  filter  hypothesis  about  jamming/spoofing  status  is  correct  (19).  The  conditional 
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probabilities  Pi (<i)  (see  Figure  1.3)  for  k  =  are  determined  by: 


PkiU)  = 


/i(t.)|a,Z(ti-i)(z.  I  at,Z,-i)Pi(t.-i) 


(2.26) 


where,  for  this  development,  the  parameters  may  represent  measurement  bias  magni¬ 
tudes  and/or  measurement  noise  strengths. 

The  numerator  of  this  equation  is  the  product  of  two  terms.  is  merely  the 

most  recent  prior  value  of  the  conditional  probability  for  the  elemental  filter;  Equation 
(2.26)  is  an  iteration.  All  K  numerator  terms  for  the  K  elemental  filters  at  tj_i  must,  of 
course,  be  available  before  the  denominator  and  consequently  the  A:*'^  probability  at  the 
current  time,  U,  can  be  calculated.  The  first  numerator  term  is  the  probability  density 
function  of  the  current  measurements  conditioned  on  both  the  assumed  parameter  values 
and  observed  past  measurements.  This  probability  density  function  is  computed  by: 


{•} 


(27r)^|A*(t,)|^®'‘^^‘^ 

{-5rn<i)Ai*(t<)rt(t<)} 


(2.27) 


where  m  is  the  measurement  vector  dimension,  and  the  filter  residual  is  given  by: 


r*(ti)  =  [z(<<)  -  Ht(t,)x)t(t,.  )]  (2.28) 

and  where  the  residual  covariance  is  computed  by  the  Ar**  elemental  filter  as: 

A,{ti)  =  [H,{U)F,itr)UliU)  +  R*(t,)]  (2.29) 

The  A:*'*  filter  residual,  rt(tj),  is  dependent  upon  the  current  measurement,  z(tj),  the  mea¬ 
surement  matrix  Ht(<j)  and  the  state  estimate,  Xt(t,“),  before  the  measurement.  The 
A;*'*  elemental  filter’s  state  estimation  error  covariance  matrix  before  the  V'  measurement, 
),  the  measurement  matrix  and  the  observation  noise  covariance  matrix,  Rt(t,),  are 
used  to  construct  Ak{ti)  (19:131).  The  Kalman  filter  equations  and  notation  were  discussed 
in  Section  2.2.1. 
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If  the  residuals  in  filter  k  are  in  consonance  with  the  filter-computed  At(ti),  the 
quadratic  in  Equation  (2.27)  is  about  equal  to  -  j,  whereas  if  the  residuals  are  much  big¬ 
ger  than  anticipated  due  to  the  wrong  parameter  hypothesis,  then  the  quadratic  is  a  much 
larger  negative  number,  so  Pk  decreases  substantially.  The  denominator  of  Equation  (2.26) 
represents  the  sum  of  the  numerator  terms  from  each  elemental  filter  computed  at  time, 
ti-  This  scales  the  numerator  to  ensure  that  the  sum  of  the  pt's  is  always  one.  However,  a 
difficulty  arises  if  the  conditional  probability  of  a  state  estimate  were  to  become  zero.  In 
such  case,  the  conditional  probability,  pt,  becomes  “locked”  at  a  zero  value  due  to  the  iter¬ 
ative  form  of  Equation  (2.27).  Real  world  changes  will  not  unlock  p*  from  zero  even  if  the 
elemental  filter  is  producing  good  state  estimates.  Therefore,  pt,  in  Equation  (2.26)  would 
remain  zero  for  all  time  because  Pfc(t,_i),  the  most  recent  value  of  the  conditional  proba¬ 
bility,  multiplies  the  other  terms  in  Equation  (2.26).  Therefore,  a  lower  bound  (threshold) 
on  Pk  should  be  set  to  preclude  such  a  zero  lock-in  conditioii , 

The  state  estimates,  x^,  are  then  scaled  (multiplied)  by  the  appropriate  weighting 
factor,  Pk,  for  each  elemental  filter  and  these  weighted  state  estimates  are  summed  for  all 
K  filters.  Summing  these  products  results  in  the  blended  or  adaptive  state  estimates: 

=  (2.30) 

fc=i 

The  covariance,  P(t/’),  of  the  blended  solution  is  estimated  by: 

(2-31) 

i=l 

The  parameter  estimates,  a*  are  calculated  in  the  same  manner  as  the  state  estimates, 
merely  by  scaling  the  assumed  value  of  the  parameter  from  each  elemental  filter  by  the 
corresponding  weighting  factor,  p*,  and  summing  these  weighted  estimates  for  all  K  filters 
according  to  the  relationship: 

K 

a(/<)  =  X;at.pt(<<)  (2.32) 

fc=i 

This  Bayesian  form  of  adaptive  estimation  is  depicted  in  Figure  1.3.  Other  approaches  are 
possible  (19)  but  will  not  be  considered  in  this  research. 
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2.4  Failure  Detection 

This  section  presents  the  theory  behind  MMAE  for  the  purpose  of  failure  detection. 
Given  a  single  KF  as  developed  in  Section  2.2,  GLR  and  chi-square  algorithms  could  be 
used  to  observe  changes  in  the  residuals.  For  significant  changes  in  the  residuals,  the  GLR 
or  chi-square  values  would  exceed  some  designated  threshold  value.  If  the  GLR  or  chi- 
square  values  exceed  that  threshold,  a  failure  is  indicated.  GLR  and  chi-square  methods 
are  mentioned  only  for  comparison  purposes  and  because  they  were  the  methods  used  by 
Vasquez  (36)  and  several  others  for  their  research.  After  completion  of  this  Section  2.4, 
only  FDI  using  MMAE  will  be  discussed  in  any  detail  in  the  remainder  of  this  thesis. 

2.4-1  FDI  via  MMAE.  This  discussion  is  presented  to  supplement  that  which  was 
presented  in  Section  2.3.  A  simplifying  fact  is  that  failure  detection  and  isolation  using 
MMAE  becomes  virtually  automatic  once  the  MMAE  has  been  developed,  as  shown  in 
Section  2.3.  Failure  detection  and  isolation  are  accomplished  merely  by  observing  the 
conditional  probabilities,  p*,  for  the  elemental  filters.  A  failure  has  been  detected  and 
identified  when  the  pk  for  the  baseline  (i.e.,  no-fail)  elemental  model  decreases  and  pk  for 
another  elemental  model  increases.  The  elemental  filter  with  the  highest  pk  contains  the 
model  which  most  closely  matches  the  truth  model  (real  world).  If  one  of  the  elemental 
filters  contained  an  assumed  parameter  value  which  exactly  matched  the  actual  failure 
magnitude,  the  failure  modeled  in  that  elemental  filter  has  therefore  been  isolated  as  the 
failure  which  occurred.  However,  if  none  of  the  elemental  filters  exactly  model  the  true 
error,  which  would  be  the  most  common  occurrence,  the  parameter  estimation  ability  could 
then  be  used  to  identify  the  failure  magnitude  which  lies  between  discrete  values  of  the 
discretized  parameter  space.  This  is  because  the  parameter  estimation  capability  would 
use  the  K  probabilities  to  weight  the  assumed  parameter  values  and  estimate  a  parameter 
value  which  would  approximate  the  true  failure  magnitude. 

2.4.2  Threshold  Selection  and  Filter  Tuning.  Revisiting  GLR  and  chi-square  FDI 
systems,  the  thresholds  must  be  selected  such  that  they  are  low  enough  to  permit  failure 
identification  and  to  minimize  delays  in  detecting  failures.  Also,  they  must  be  low  enough 
to  prevent  missed  alarms  caused  by  GLR  or  x  values  dropping  below  the  threshold  while 
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an  actual  failure  is  still  present.  Finally,  the  thresholds  must  be  high  enough  to  prevent 
false  alarms  caused  by  variations  in  the  GLR  and  x  values.  These  variations  may  result 
from  aircraft  maneuvers  or  unpredictable  changes  in  the  random  measurement  noise  or 
from  other  unmodelled  effects.  If  the  filter  is  tuned  sufficiently  for  both  good  tracking 
and  enhanced  failure  detection,  such  variations  due  to  aircraft  maneuvers  or  unpredictable 
changes  in  the  random  measurement  noise  define  a  relatively  low  bound  set  by  the  noise 
magnitude.  Analogously,  assuming  adequate  tuning  of  the  elemental  filters,  thresholds  for 
FDI  with  MMAE  are  accomplished  by  observing  the  elemental  filter  probabilities  ranging 
from  0  to  1  such  that  failures  are  detected  and  identified  by  the  elemental  filter  probabilities. 
Lower  bound  on  the  probabilities  must  be  set  appropriately,  high  enough  to  prevent  zero 
lock-in  but  low  enough  still  to  allow  meaningful  probability  calcultions.  If  the  lower  bound 
is  too  high,  probabilities  below  that  value  are  impossible  to  attain,  and  the  adaptation  is 
incapacitated  because  too  little  of  the  total  probability  is  still  available  to  be  moved  around 
among  the  elemental  filters. 

When  tuning  a  Kalman  filter,  a  major  tradeoff  must  be  made  to  meet  the  goals  for 
state  estimation  and  failure  detection.  Adjusting  the  process  noise  and  measurement  noise 
values  to  enhance  FDI  capabilities  may  degrade  state  estimation  and  vice  versa.  This  can 
be  seen  in  Equation  (2.29)  by  realizing  that  reductions  in  R  will  cause  reductions  in  A. 
The  consequence  is  that  reductions  in  R  will  cause  an  increase  in  the  ratio  of  Q  eigenvalues 
to  R  eigenvalues,  resulting  in  a  conservatively  tuned  filter  which  may  yield  performance 
degraded  from  the  performance  achievable  when  the  filter  is  tuned  specifically  for  best 
state  estimation  performance.  Decreasing  R  and  increasing  Q  could  yield  the  same  tuning 
result,  though,  since  the  steady  state  gain,  K,  depends  on  the  ratio  of  Q  and  R  eigenvalues. 
Observation  of  the  state  estimates,  measurement  residuals  and  the  residual  covariance  A(t,) 
will  indicate  a  point  of  diminishing  returns  using  this  technique.  The  elemental  filters  can 
be  tuned  either  via  multi-run  MSOFE  simulations  or  tuned  using  MMSOFE  with  different 
tuning  values  used  in  the  different  elemental  filters. 
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Figure  2.1.  MMAE-Based  Control 


2.5  Stochastic  Adaptive  Control 

In  a  generic  sense  of  MMAE  simulations,  MSOFE  could  be  used  to  run  each  elemental 
filter  independently  from  initial-time  to  final-time,  with  the  MMAE  aspects  (state  and 
covariance  blending,  parameter  estimation,  and  etc.)  being  conducted  post-process.  This 
is  not,  however,  really  feasible  if  the  MMAE-optinaized  values  are  used  in  calculations 
which  affect  the  elemental  filter  calculations  during  the  ensuing  propagation  and  update 
cycles.  With  that  in  mind  and  although  not  necessarly  an  exhaustive  list,  several  types  of 
stochastic  adaptive  control  which  are  discussed  in  Chapter  15  of  Maybeck’s  Volume  3  (20) 
are  of  interest  with  regards  to  MMSOFE. 

Closed-loop  forms  of  stochastic  adaptive  control  which  use  MMAE  to  estimate  states 
and  parameter  values  optimally  require  a  true  MMAE  structure  as  available  with  MM¬ 
SOFE.  These  parameter  values  are  then  used  in  an  adaptive  control  law  as  shown  in  Figure 
2.1,  known  as  an  MMAE-Based  Controller.  In  using  this  type  of  controller,  Gc[ti;a(t,^)] 
is  determined  as  based  on  the  parameter  estimate  a  which  is  assumed  to  remain  con¬ 
stant  until  the  next  measurement  update.  Finally,  Maybeck  (20)  discusses  Multiple  Model 
Adaptive  Controllers  (MMAC)  as  a  type  of  stochastic  adaptive  controller  which  consists 
of  a  bank  of  Linear  Quadratic  Gaussian  (LQG)  controllers  in  place  of  just  the  elemental 
filter  of  MMAE  (see  Figure  2.2).  Again  with  MMAC,  the  MMAE  structure  is  used  to 
calculate  probabilities  for  use  as  weighting  factors  to  blend  the  control  signals,  u^,  from 
the  elemental  controUers  optimally. 


Figure  2.2.  Multiple  Model  Adaptive  Controller 


2.6  Distributed  Kalman  Filtering 

Distributed  Kalman  filter  (DKF)  simulations  represent  another  problem-type  which 
could  well  be  accomplished  with  a  slightly  modified  MMSOFE.  Operationally,  DKF  would 
run  its  various  KFs  on  separate  processors  with  a  master  (central)  filter  blending  the  out¬ 
put  of  the  individual  filters.  However,  for  simulation  purposes,  MMSOFE  as  conceived 
for  MMAE  in  this  thesis  could  be  modified  only  slightly  to  implement  a  DKF  simula¬ 
tion.  In  such  a  configuration,  even  raaster-filter-calculated  values  could  be  fed  back  to 
the  distributed  filters.  One  type  of  DKF  implementation  is  depicted  in  Figure  2.3.  In  this 
DKF  each  local  (elemental  in  MMAE)  filter  can  have  a  unique  model  and  receive  measure¬ 
ments  from  different  local  sensors  (the  same  measurements  feed  into  all  elemental  filters 
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Matter  Feedback 


dx 


Figure  2.3.  One  Possible  Distributed  Kalman  Filter  Estimator 


in  MMAE).  These  local  filters  calculate  the  estimates  of  their  unique  set  of  states  with 
related  covariances,  both  of  which  are  forwarded  from  aU  of  the  local  filters  to  the  master 
filter.  The  master  filter’s  function  is  to  combine  or  integrate  all  of  these  state  estimates  via 
some  problem-specific  combinatorial  algorithm  or  even  via  a  Kalman  filter.  This  master 
filter  may  generate  some  single  master  feedback  or  set  of  master  feedbacks  which  are  fed 
back  to  the  reference  system  and/or  local  filters. 


2. 7  Summary 


This  chapter  has  provided  the  theoretical  basis  for  the  remaining  chapters.  MMAE, 


FDI  techniques,  MMAE-based  control  and  MMAC,  and  distributed  Kalman  filtering  have 
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been  discussed  in  Chapter  II.  Chapter  III  addresses  the  actual  software  implementation  of 
MMA£  which  was  developed  in  this  thesis. 
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III.  MMSOFE  Development 


5.1  Oveniew 

This  chapter  discusses  the  aspects  of  MSOFE  which  are  germane  to  MMSOFE  devel¬ 
opment,  including  a  brief  description  of  the  satellite  orbit  problem.  A  cursory  discussion 
of  MSOFE  and  a  more  detailed  description  of  the  MMSOFE  adaptation  of  MSOFE  are 
made,  including  the  actual  modifications  made  to  MSOFE.  The  validation  of  MMSOFE  is 
discussed  next,  and  finally,  MMSOFE  capabilities. 

5.2  Satellite  Orbit  Problem 

The  satellite  orbit  problem  is  fully  described  in  Maybeck’s  Volumes  I  and  II  (18,  19) 
and  in  the  User's  Manual  for  MSOFE  (24),  but  since  it  is  the  test  problem  for  MMSOFE 
development,  it  warrants  some  minimal  discussion  in  this  thesis.  The  problem  is  con¬ 
structed  for  MSOFE  such  that  the  simulated  satellite  is  already  in  an  orbit  with  an  orbit 
radius  of  one  “radius  unit”  from  the  center  of  the  earth.  “Radius  units”  and  “time  units” 
are  used  as  non-dimensional  units  merely  as  a  simplifying  measure.  The  MSOFE  user  can 
then  implement  a  Hohmann  orbit  transfer  as  well  as  modify  the  initial  conditions,  noise 
variances,  and  choose  MSOFE  modes  and  parameters  via  the  input  data  file  MSOFEJN. 
The  user  can  specify  the  times  for  the  Hohmann  transfer  actions  and  the  ratio  of  the  orbit 
radii  before  and  after  completion  of  the  Hohmann  transfer. 

As  mentioned  in  Section  3.3,  the  failures  are  added  into  the  problem  in  the  sub¬ 
routine  HRZSYS.  At  the  times  specified,  the  failures  were  induced  in  HRZSYS  with  the 
measurement  bias/ramp  being  added  directly  to  the  truth  measurement  and  the  mea¬ 
surement  noise  variance  bias/ramp  added  to  the  baseline  (no-failure)  measurement  noise 
variance.  Additionally,  the  Hohmann  transfer  was  not  implemented  for  these  simulations 
during  software  development,  in  order  to  accentuate  the  performance  of  the  independent 
elemental  filters,  the  performance  of  the  blended  estimation,  and  that  of  the  probability 
and  parameter  calculations. 
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3.S  MSOFE  Description 


As  mentioned  in  Chapter  I,  MSOFE  provides  a  tool  for  both  Monte  Carlo  simulations 
and  covariance  analysis  of  linear  or  extended  Kalman  filters.  The  version  of  MSOFE  used 
for  this  development  was  Version  1.5,  which  was  current  as  of  June  1991.  The  MSOFE 
FORTRAN  code  is  composed  of  two  source  code  files,  MSOFE.f  and  USOFE.f.  MSOFE.f, 
containing  59  subroutines,  comprises  the  “core-code”  for  MSOFE  and  users  should  rarely 
need  to  access  or  modify  MSOFE.f.  USOFE.f,  containing  17  subroutines,  is  the  user- 
interface  source  code.  While  some  of  these  USOFE.f  subroutines  remain  as  part  of  the 
MSOFE  code  even  when  they  aren’t  being  used  for  a  ^ven  application  and  are  just  empty 
shells,  they  are  retained  so  as  to  ease  the  adaptation  to  other  specific  filtering  problems  in 
which  they  are  required. 

In  addition  to  “adapting”  USOFE.f  for  the  specific  problem,  the  input  and  output 
hies  which  are  the  other  user-interfaces  must  also  be  “adapted.”  MSOFE  input  and  output 
files  are  shown  in  Table  3.1. 

Table  3.1.  MSOFE  Input/Output  Files 


File 

Name 

FORTRAN 
Unit  # 

Description 

INPUT  FILES 

MSOFEJN 

5 

User  controls  & 
initial  conditions 

FLIGHT 

3 

Externally  produced  trajectory 

OLDEND 

9 

Data  from  last  run 

KFGAIN 

11 

Kalman  gain  series 

OUTPUT  FILES 

CTOM 

2 

Continuous-time  data 

DTOM 

4 

Discrete-time  data 

MSOFE-OUT 

6 

Text  file 

10 

Data  for  next  run 

KFGAIN 

12 

Kalman  gain  series 

ERRS 

8 

Error  messages 

INPUT/OUTPUT 

scratch 

7 

Scratch  space 

The  FORTRAN  unit  numbers  are  merely  numbers  assigned  to  each  particular  file  when 
the  file  is  opened,  by  which  FORTRAN  refers  to  files  during  execution.  To  explain  the 
files  referenced  in  Table  3.1  further: 

MSOFE  JN :  Contains  controls  to  direct  MSOFE  operation  and  parameters  to  define  the 
problem  and  specify  initial  conditions.  MSOFE  JN  is  broken  into  five  groups:  Prob¬ 
lem  Title,  Control  Parameters,  Printer  Plots  Specifications,  Initial  Conditions  and 
User  Data. 

CTOM  &  DTOM:  Serves  as  primary  data  output  files  for  use  by  plotting  or  other  post¬ 
process  software. 

FLIGHT:  Provides  the  input  trajectory  data  for  many  MSOFE  problems,  such  as  for 
Vasquez’s  (36)  research  in  which  the  “Profile  Generator”  (PROFGEN)  software  was 
used  to  generate  the  FLIGHT  trajectory  file.  FLIGHT  could  alternatively  be  created 
via  other  software  available  to  the  user.  However,  the  satellite  orbit  problem  used 
in  MMSOFE  development  has  its  trajectory  generation  software  built  into  USOFE.f 
code. 

KFGAIN:  Can  be  generated  by  MSOFE  and  then  used  as  input  for  future  MSOFE  runs 
to  specify  the  precomputed  Kalman  gain  vector  during  the  simulation. 

MSOFE-OUT:  Generated  by  MSOFE.f,  provides  an  ASCII  record  of  control  values,  initial 
conditions,  and  any  other  problem-specific  data  provided  in  MSOFEJN.  Could  also 
contain  additional  problem-specific  data  as  specified  by  the  user  in  USOFE.f.  Also 
contains  error  messages  and  printer  plots  if  requested. 

NEWEND  &  OLDEND:  Final  condition  data  is  saved  to  NEWEND  at  the  end  of  an 
MSOFE  run  for  use  as  initial  conditions  in  OLDEND  for  subsequent  MSOFE  simu¬ 
lations.  The  change  of  name  from  NEWEND  to  OLDEND  is  handled  externally  by 
the  user. 

ERRS:  Contains  any  MSOFE-generated  error  messages  produced  during  or  causing  the 
termination  of  a  run. 
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Plots  of  the  data  in  CTOM  and  DTOM  can  be  made  using  any  of  a  variety  of 
plotting  packages.  In  this  version  of  MSOFE,  CTOM  and  DTOM  were  generated  in  a 
format  compatible  as  input  for  MPLOT  which  is  described  in  the  MSOFE  user's  manual 
(24).  MPLOT  was  then  used  to  create  MatriXx-compatible  data  files  for  plotting  with 
MatriXx  (10). 

Because  a  thorough  discussion  of  MSOFE  is  available  in  the  User’s  Manual  for 
MSOFE  (24),  that  effort  will  not  be  repeated  here,  but  rather  only  those  parts  of  MSOFE 
which  were  appreciably  altered  or  which  bear  special  significance  in  MMSOFE  develop¬ 
ment  will  be  discussed  in  any  depth.  However,  to  clarify  the  modifications  to  MSOFE, 
some  flowchart-like  diagrams  are  included  in  this  chapter.  Figures  3.1  -  3.7.  First,  the 
MSOFE  top-level  executive  program,  referred  to  as  MAIN,  calls  several  subroutines  as 
shown  in  Figure  3.1. 

These  flowchart-like  diagrams.  Figures  3.1  -  3.7,  are  obviously  not  conventional 
flowcharts.  The  MMSOFE  diagrams,  Figures  3.4  to  3.7,  were  made  first  and  the  MSOFE 
diagrams,  Figures  3.1  to  3.3,  were  made  by  deleting  the  non-applicable  parts.  This  results 
in  some  of  the  MSOFE  diagrams  appearing  slightly  sparse  but  allows  for  better  compar¬ 
isons  between  the  MSOFE  and  MMSOFE  diagrams.  The  “Note”  in  these  diagrams  also 
requires  some  additional  explanation.  These  figures  are  not  intended  as  detailed  flowcharts 
of  program  code.  Rather,  they  show  the  flow  and  order  of  subroutine  calls  as  they  appear 
in  MSOFE  or  MMSOFE  code  and  the  related  branch  statements.  Except  for  a  few  loops 
included  for  clarity’s  sake,  most  of  the  code  is  omitted  from  these  MSOFE  and  MMSOFE 
diagrams  (Figures  3.1  -  3.7).  Recalling  that  the  “M”  in  MSOFE  stands  for  Multimode, 
the  “If _ ”  statements  contained  in  some  function  blocks  or  subroutine  call-ellipses  in¬ 
dicate  that  action  or  call  statement  is  conditional  on  the  “If _ ”  statement  being  true. 

Some  of  these  “If _ ”  conditions  vary  between  passes  through  the  code  during  a  single 

Monte  Carlo  “run,”  or  for  that  matter,  for  different  “runs”  of  a  Monte  Carlo  simulation. 
For  excimple,  if  MSOFE  is  run  in  “dual”  mode,  both  covariance  analysis  and  Monte  Carlo 
simulation,  the  code  for  both  is  implemented  during  the  first  Monte  Carlo  “run,”  but  the 
code  specific  only  to  the  covariance  analysis  isn’t  implemented  during  subsequent  Monte 
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Figure  3.1.  MSOFE  MAIN  Program 

Carlo  “runs.”  Furthermore,  since  the  MSOFE  User’s  Manual  (24)  describes  the  function 
of  all  of  the  MSOFE  subroutines,  they  won’t  be  described  in  this  thesis. 

The  blocks  in  Figures  3.1  and  3.4  containing  lists  of  “Subroutines  Called”  mean  that 
these  subroutines  are  called  by  the  subroutine  that  is  named  in  the  ellipse  directly  above 
the  box.  The  major  functions  of  these  subroutines  are  as  follows: 

•  BLKDAT:  “Block  Data”  -  Initializes  all  common  areas  with  “DATA”  statements. 


3-5 


•  INITPR:  “Initialize  Problem”  -  Called  at  the  beginning  of  the  simulation  to  read 
in  problem-specific  data  from  MSOFE  JN  and  to  initialize  all  variables  common  to 
all  Monte  Carlo  runs. 

-  MACHIN:  “Machine”  -  Computes  several  computer-specific  constants  which 
characterize  the  computer. 

-  RD12:  “Read  Groups  1  &  2”  -  Reads  Group  1  and  Group  2  data  from 
MSOFEJN  and  writes  header  information  to  MSOFE.OUT. 

-  VALPAR:  “Validate  Parameters”  -  Examines  parameter  from  SIZES  to  ensure 
they  fall  within  range  restrictions. 

-  TABLEl:  “Table  1”  -  Writes  TABLE  1  which  contains  a  list  of  parameters 
from  the  SIZES  file  and  MSOFEJN  Group  2  data  to  MSOFE.OUT. 

—  TABLE2:  “Table  2”  -  Writes  computer-specific  data  to  MSOFE.OUT. 

-  VALDAT:  “Validate  Data”  -  Examines  Group  2  parameters  to  ensure  they  fall 
within  range  restrictions  and  to  check  filter  and  composite  covariance  matrices 
for  positive  definiteness. 

-  04PLOT:  “Output  to  Unit  #2  &  #4”  -  Writes  data  to  CTOM  and  DTOM  for 
postrun  processing. 

-  07PLOT:  “Output  for  Plotting”  -  Executive  routine  for  making  printer  plots 
of  user-specified  variables. 

-  012KF:  “KF  Output  to  Unit  #12”  -  Writes  the  Kalman  gain  to  the  file  KF- 
GAIN. 

•  INITRN:  “Initialize  Run”  -  Called  at  the  beginning  of  each  simulation  run  to  read 
in  run-specific  data  and  to  initialize  run-specific  variables. 

-  GETP:  “Get  P  (covariance)”  -  Reads  MSOFEJN  to  acquire  the  initial  values 
for  the  filter,  filter  error,  and  truth  covariance  matrices. 

—  RDHEDR:  “Read  Header”  -  Reads  the  header  of  the  FLIGHT  trajectory  file. 
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—  SPAN:  “Span  T”  -  Reads  the  data  from  FLIGHT  trajectory  file  and  manipu¬ 
lates  data  so  the  simulation  times  are  synchronized  with  trajectory  data. 

-  INTERP:  “Interpolate”  -  Interpolates  trajectory  data  to  the  time  series  being 
used. 

-  MFIL:  “Matrix  Fill”  -  Fills  a  matrix  with  a  scalar  values. 

-  VMOVE:  “Vector  Move”  -  Copies  one  vector  into  another. 

-  UCOREL:  “Correlated  Gaussian  Sample  Vector”  -  Constructs  a  Gaussian  sam¬ 
ple  vector  with  zero  mean. 

-  TRJSYS:  “Trajectory,  System”  -  Calculates  user-specified  trajectory  data. 

-  FSMAP:  “Filter /System  Mapping”  -  Computes  the  matrix  AFS  which  maps 
truth  model  states  into  equivalent  filter  states. 

-  PCINIT:  “Pc  Initialize”  -  Initializes  the  composite  covariance  matrix. 

-  PFINIT:  “Pf  Initialize”  -  Initializes  the  filter  covariance  matrix. 

-  USRIN:  “User  Input”  -  Reads  data  as  specified  by  the  user. 

-  SNOYS:  “System  Noise”  -  Provides  random  noise  to  the  truth  state  vector. 

-  VALDAT :  See  description  under  INITPR  above. 

—  DERIV:  “Derivative”  -  Called  to  initialize  the  derivatives  of  the  truth  and  filter 
states,  as  well  as  of  the  filter  and  composite  covariances. 

-  OUT:  “Output”  -  Controls  the  outputs  which  are  specified. 

•  SIMRUN;  “Simulation  Run”  -  The  top-level  controller  for  the  time  cycles  of  simu¬ 
lation  runs. 

-  RDGAIN:  “Read  Gain”  -  Reads  the  next  measurement  time  and  the  Kalman 
gain  from  the  KFGAIN  file. 

-  TNEXTM:  “Time  of  Next  Measurement”  -  Calculates  the  next  measurement 
time  for  each  measurement  type. 

—  TNEXTE:  “Time  of  Next  Event”  -  Determines  the  time  of  the  next  event. 
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-  ADVANS:  “Advance”  -  ContrHs  the  propagation  between  events  via  numerical 
integration. 

-  SPAN:  See  description  under  INITRN  above. 

—  SNOYS:  See  description  under  INITRN  above. 

-  KUTMER:  “Kutta  Merson”  -  Integrates  a  set  of  first-order  ordinary  linear  or 
non-linear  differential  equations. 

—  DERJV:  See  description  under  INITRN  above. 

•  OUT:  See  description  under  INITRN  above. 

•  MEASUP:  “Measurement  Update”  -  Performs  the  measurement  update  for  each 
scalar  measurement  one  at  a  time. 

-  HRZSYS:  “System  H,  R,  Z”  -  Supplies  the  truth  measurement  (Z),  the  truth 
measurement  noise  variance  (R),  and  the  partial  derivatives  of  the  true  mea¬ 
surement  with  respect  to  the  filter  states  and  the  system  states  (H), 

-  HRZFIL:  “Filter  H,  R,  Z”  -  Supplies  the  filter-predicted  measurement  (Z),  the 
filter  measurement  noise  variance  (R),  and  the  partial  derivatives  of  the  filter 
measurement  with  respect  to  the  filter  states  and  the  system  states  (H). 

-  PSQRT:  “Square  Root  of  the  System  Covariance”  -  Calculates  the  square  root 
of  the  covariance  matrix  in  terms  of  the  unitary  (U)  and  diagonal  (D)  matrices. 

—  UDMESC:  “Composite  UD  Measurement”  -  Performs  the  measurement  update 
of  the  composite  covariance  U  and  D  matrices. 

-  UDMESF:  “Filter  UD  Measurement”  -  Performs  the  measurement  update  of 
the  filter  covariance  U  and  D  matrices. 

—  UDSQC;  “Composite  UD  Squared”  -  Forms  the  composite  covariance  matrix 
in  upper  triangular  form  from  the  U  and  D  factors. 

—  UDSQF:  “Filter  UD  Squared”  -  Forms  the  filter  covariance  matrix  in  upper 
triangular  form  from  the  U  and  D  factors. 
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•  FDBACK;  “Feedback”  -  Updates  the  filter  and  system  states  and/or  the  composite 

covariance  matrix  with  regard  to  the  discrete  feedback  correction. 

-  XFDBAK:  “State  Feedback”  -  Updates  truth  and  filter  states  with  discrete 
feedback  correction. 

-  DFDBAK;  “D  Feedback”  -  Determines  the  matrices  used  to  pre-multiply  the 
truth  or  filter  state  vectors  to  provide  discrete  corrections  to  the  states. 

BLKDAT  (block  data)  contains  the  default  initialization  DATA  statements  to  initial¬ 
ize  the  values  of  variables,  in  contrast  to  PARAMETER  statements  which  set  the  values  of 
parameters.  The  BLKDAT  “call”  is  handled  as  “EXTERNAL  BLKDAT”,  not  as  a  normal 
subroutine  call.  Parameter  values  cannot  be  changed  in  the  program  but,  of  course,  vari¬ 
ables  can  be  changed.  Conversely,  parameters  can  be  used  to  dimensionalize  vectors  and 
matrices  in  declaration  statements  but  variables  cannot.  This  “EXTERNAL  BLKDAT” 
statement  has  the  effect  of  placing  all  the  DATA  statements  contained  in  BLKDAT  right 
there  where  the  EXTERNAL  statement  is  found  without  having  all  the  DATA  statements 
actually  present  in  the  top-level  executive  routine. 

The  MSOFE  subroutines  INITPR,  INITRN,  SIMRUN,  and  OUT  are  all  called  from 
MAIN  and  are  all  germane  to  MMSOFE  development.  INITPR  is  the  subroutine  which 
initializes  the  variables  and  flags  for  the  entire  MSOFE  simulation  and  is  only  executed 
once  for  MSOFE.  INITRN  initializes  some  variables  and  flags  for  each  run  in  a  Monte 
Carlo  MSOFE  simulation  and  is  executed  at  the  beginning  of  each  run.  SIMRUN,  the 
MSOFE  second-level  executive  subroutine,  calls  a  number  of  other  subroutines,  as  shown 
in  Figure  3.2.  SIMRUN  is  responsible  for  timing  and  scheduling  propagations,  measure¬ 
ment  updates,  feedback,  and  judiciously  calling  OUT  with  the  appropriate  argument.  The 
OUT  subroutine  controls  MSOFE  scheduled  output,  special  prints,  and  a  few  other  special 
functions.  The  OUT  integer  argument,  lOUT,  shown  in  the  diagrams  with  the  OUT  calls, 
specifies  the  event  which  has  occurred  and  determines  OUT’s  internal  flow. 

With  respect  to  MMSOFE  development,  the  call  to  MEASUP  (Figure  3.3)  is  the 
most  notable  subroutine  call  in  SIMRUN.  It  is  the  most  significant  because,  in  MMSOFE, 
MEASUP  is  where  nearly  all  of  the  calls  to  the  MMSOFE  MMAE  subroutine  are  made 
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(see  Figures  3.4  -  3.7).  The  MMSOFE  MMAE  subroutine,  the  code  of  which  is  presented 
in  Figure  A.l  and  which  is  discussed  in  Section  3.4.2.3,  performs  nearly  all  of  the  multiple 
model  related  calculation  and  also  writes  to  the  output  file,  BLCNTOM  (data  file  con¬ 
taining  the  blended  or  MMAE-related  output).  So,  the  filename  BLCNTOM  was  derived 
from  CTOM,  the  MSOFE  output  file,  with  BL  =  blended  and  N  merely  refering  to  N  files 
for  N  elemental  filters. 

Subroutine  SIMRUN  handles  the  time  progression  of  measurements  and  other  events 
for  MSOFE  (TNEXTM,  TNEXTE)  and  the  subroutine  calls  for  propagation  (ADVANS), 
measiirement  update  (MEASUP),  subroutine  calls  for  the  application  of  discrete  feedback 
(FDBACK),  and  the  outputs  at  appropriately  scheduled  times  (OUT).  The  subroutine 
SIMRUN  performs  all  of  these  functions  for  all  Monte  Carlo  runs  before  returning  control 
to  the  MSOFE  MAIN  program. 

The  MEASUP  subroutine,  which  is  called  by  SIMRUN,  controls  the  calculations 
related  to  the  measurement  update  by  calling  subroutines  HRZSYS  and  HRZFIL.  In 
MSOFE’s  subroutine  MEASUP,  HRZSYS  and  HRZFIL  are  the  subroutine  calls  which 
are  most  pertinent  to  MMSOFE  development.  These  two  subroutines  reside  in  USOFE.f 
and  are  therefore  intended  for  user  interface.  They  build  the  partial  derivative  matrices 
which  are:  derivative  of  the  true  measurement  with  respect  to  the  truth  states  (HSYSS), 
derivative  of  the  true  measurement  with  respect  to  the  filter  states  (HSYSF),  derivative 
of  the  filter  measurement  with  respect  to  the  truth  states  (HFILS),  and  derivative  of  the 
filter  measurement  with  respect  to  the  filter  states  (HFILS).  These  are  used  in  MEASUP 
to  build  the  measurement  observation  vectors.  HRZSYS  and  HRZFIL  also  build  the  truth 
and  filter  measurements  respectively,  including  measurement  noise  addition  for  the  truth 
measurement.  For  the  orbit  estimation  problem  as  used  for  this  thesis,  HRZSYS  is  where 
the  failures  given  in  Table  3.1  were  instilled  into  the  true  range  measurement  or  true  mea¬ 
surement  noise  variance.  Likewise,  HRZFIL  was  used  to  add  the  elemental  filter-assumed 
measurement  signal  biases  or  measurement  noise  strength  alterations  for  both  MSOFE 
and  MMSOFE.  Subroutine  MEASUP  also  controls  the  calls  to  subroutines  UDMESC  and 
USMESF  which  perform  the  U  —  D  covariance  factorization  calculations  for  updating  the 
covariances  for  the  measurement  updates,  and  it  controls  the  calls  to  the  UDSQC  and 
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Figure  3.3.  MSOFE  MEASUP  Subroutine 
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3.4  MMSOFE  Adaptation  from  MSOFE 

UDSQF  subroutines  which  form  the  appropriate  covariance  matrice  elements  and  stores 
them  in  a  vector  containing  the  upper  triangular  matrix  elements. 

After  it  had  been  determined  that  SIMRUN  would  be  the  top-level  executive  subrou¬ 
tine  for  MMSOFE,  the  methodology  followed  for  MMSOFE  development  was  as  discussed 
in  Section  1.6.4.  First,  The  MMSOFE  input  and  output  files  will  be  described.  Then, 
the  new  subroutines  developed  for  MMSOFE  (EMFIL,  INVERT,  MMAE,  and  OPFILE) 
are  discussed,  followed  by  discussions  of  the  MMSOFE  versions  of  subroutines  modified 
from  MSOFE  subroutines  (MAIN,  SIMRUN,  amd  MEASUP)  and  the  MMATH  subroutine 
which  was  borrowed  and  modified  from  Version  1.5  of  MSOFE.f,  current  in  February  1993 
and  named  CSOFE.f.  CSOFE.f  is  actually  the  version  of  MSOFE  described  in  the  User’s 
Manual  for  MSOFE  (24).  Finally  the  “vectorization”  of  MSOFE  scalar  variables  (mak¬ 
ing  elemental  filter-specific  or  blended  output-specific  variables)  and  “matricization”  of 
MSOFE  vectors  will  be  discussed,  but  both  will  be  referred  to  as  “vectorization”. 

3,4.1  Input/Output  Files.  Because  more  input  and  output  files  are  required  for 
MMSOFE  than  for  MSOFE  (see  Table  3.1),  a  new  subroutine,  OPFILE  (Section  3.4.2.4  and 
Figure  A.2),  was  devised  to  accommodate  automatically  opening  and  appropriately  naming 
the  input  and  output  files.  OPFILE  assigns  names  and  FORTRAN  unit  file  numbers  to 
the  MMSOFE  input  and  output  files  as  required  by  the  number  of  elemental  filters  being 
utilized.  The  names  and  unit  numbers  of  the  hies  opened  by  OPFILE  correspond  somewhat 
to  the  files  used  by  conventional  MSOFE.  The  following  table.  Table  3.2,  shows  that 
correspondence  with  MSOFE  and  indicates  the  file  naming  convention  used  by  MMSOFE. 
There  are  several  points  of  explanation  necessary  with  regard  to  Table  3.2.  The  “xx” 
refers  to  the  total  number  of  elemental  filters  used.  If  the  computer  being  used  could 
handle  such  a  large  problem  quickly  enough,  OPFILE  would  allow  MMSOFE  to  process 
an  MMAE  simulation  with  up  to  98  elemental  filters  (i.e.  xx  =  98).  BLCNTOM  is  the 
output  file  for  the  MMAE  blended  outputs,  probabilities,  and  parameter  estimates,  and  all 
of  their  related  statistics.  BLDSTOM  was  not  used  for  this  thesis  effort,  but  it  is  opened 
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Table  3.2.  MMSOFE  Input/Output  Files 


MSOFE  File 
Name 


MSOFE 
File  Unit  # 


MMSOFE  File 
Names 


FLIGHT 


OLDEND 


KFGAIN 


OlINPUT  - 

501  - 

xxINPUT 

5xx 

&  BLINPUT 

Sxx-l-l 

unchanged 

3 

OLEND 

9 

unchanged 

11 

CTOM 


DTOM 


MSOFE.OUT 


NEWEND 


KFGAIN 


ERRS 


OICNTOM  - 

201  - 

xxCNTOM 

2xx 

&  BLCNTOM 

2xx4-l 

OIDSTOM  - 

401  - 

xxDSTOM 

4xx 

&  BLDSTOM 

4XX-I-1 

OlMDATA  - 

601  - 

xxMDATA 

6xx 

&  BLMDATA 

6xx-i-l 

unchanged 

10 

unchanged 

12 

OlERRRS  - 

801  - 

xxERRRS 

8xx 

&  BLERRRS 

8XX-1-1 

01 scratch  - 

701  - 

xxscratch 

7xx 

&  BLscratch 

7XX-H 
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by  MMSOFE  and  could  be  used  by  an  MMSOPE  user.  The  distinctions  “C  =  continuous” 
in  CTOM  and  the  “D  =  discrete”  in  DTOM  are  determined  by  the  header  information,  not 
by  the  OPEN  statement,  and  those  distinctions  are  maintained  in  MMSOFE  _CNTOM 
and  —OSTOM.  Similarly,  the  BLINPUT,  BLMDATA  and  BLERRRS  files  are  not  used 
but  are  also  opened  and  available  for  use.  The  BLINPUT  file  was  specifically  left  for  use 
as  an  input  file  for  data  which  might  be  required  for  some  of  the  other  applications  such 
as  MMAE-based  stochastic  based  control,  MMAC,  or  DKF. 

One  other  caviate  to  the  MMSOFE  file  structure  has  relevance  only  when  MMSOFE 
is  used  to  simulate  a  single  filter.  When  there  is  only  one  filter  instead  of  multiple  filters, 
the  “prefix”  will  be  absent  from  all  filenames  and  the  “BL”  files  aren’t  utilized  either.  For 
example,  the  output  files  _CNTOM  and  _I)STOM  would  be  just  CNTOM  and  DSTOM, 
and  there  would  be  no  BLCNTOM  nor  BLDSTOM.  Similarly,  —INPUT  would  be  just 
INPUT  and  there  would  be  no  BLINPUT.  This  feature  was  designed  to  minimize  memory 
use  when  performing  single  filter  simulations. 

The  —INPUT  files  are  identical  to  MSOFE  JN  except  for  several  parameters  added  at 
the  end  of  the  —INPUT  files  to  ease  the  modification  of  problem-specific  and/or  elemental 
filter-specific  variables.  This  is  not  to  say,  however,  that  all  of  the  values  input  via  the 
—INPUT  files  can  be  different  for  every  filter.  An  MMSOFE  user  might,  for  instance, 
want  to  vary  the  number  of  measurements  (NM)  between  different  elemental  (local)  filters 
in  a  DKF  application  or  vary  the  initial  (TMI)  or  final  times  (TMF)  of  measurement 
availability  for  some  application,  but  these  variables  which  are  read  from  the  —INPUT  files 
would  first  require  vectorization.  While  MMSOFE  could  be  further  modified  to  facilitate 
an  implementation  for  some  particular  application,  because  MMSOFE  reads  from  the 
—INPUT  files  in  numerical  order,  MMSOFE  will  retain  the  last  value  read  for  any  variable 
which  hasn’t  been  vectorized,  namely  the  value  from  the  xxINPUT  file  (highest  numbered 
—INPUT  file).  Those  variables  which  have  been  vectorized  are  discussed  in  Section  3.4.4. 
For  the  satellite  orbit  problem  test-case,  the  variables  added  at  the  end  of  the— INPUT 
files  include: 
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1.  RNERR  is  the  number  corresponding  to  the  failure  type  implemented  in  HRZSYS 
and  is  a  scalar,  remaining  constant  for  all  elemental  filters. 

2.  MINPROB  is  the  lower  bound  value  for  the  probabilities  to  prevent  the  “locked  in” 
situation  described  in  Section  2.3  and  is  also  a  scalar. 

3.  ZFBIAS  is  a  vector  of  values  containing  the  hypothesized  measurement  signal  bias 
magnitudes  corresponding  to  the  individual  elemental  filters. 

4.  RFNOIS  is  a  vector  of  values  corresponding  to  the  hypothesized  measurement  noise 
variance  magnitudes  applied  to  the  different  elemental  filters. 

The  calls  to  subroutine  OPFILE  to  generate  the  xx  files  corresponding  to  MSOFE’s 
files  FLIGHT,  OLDEND/NEWEND,  and  KFGAIN  were  included  in  MMSOFE,  but  they 
weren't  implemented  in  this  thesis  because  none  of  these  four  files  were  required  for  the 
orbit  estimation  problem.  If  an  MMSOFE  user  wanted  to  conduct  a  study  using  such  files 
specific  to  each  elemental  filter,  these  call  statements  could  be  easily  activated. 

S.4-S  New  Subroutines.  The  two  completely  new  subroutines  are  included  in  Ap¬ 
pendix  A;  MMAE  in  Figure  A.l  and  OPFILE  in  Figure  A.2. 

3.4.2,  t  EMFIL  Subroutine.  The  EMFIL  subroutine  computes  the  MMAE 
blended  state  estimation  error  vector,  E^uen<ud  given  by: 

Etiended  =  ^Blended  ~  AFS  •  \Truth  (3.1) 

where  'Kaunded  represents  the  MMAE-blended  state  estimates,  AFS  is  the  filter/system 
mapping  matrix  computed  by  subroutine  FSMAP  in  USOFE.f,  and  Xrruth  is  the  true  state 
vector.  EMFIL  was  adapted  from  the  MSOFE.f  subroutine  EFFIL  which  finds  the  state 
estimation  error  for  the  single  filter  and  for  the  elemental  filters  in  MMSOFE. 

3.4.2. 2  INVERT  Subroutine.  The  INVERT  subroutine  takes  the  inverse  of  a 
square  matrix  of  any  dimension  using  the  method  of  Gaussian  elimination.  It  was  written 
in  1988  by  students  at  the  Air  Force  Institute  of  Technology  advised  by  Dr.  Peter  Maybeck 
on  projects  sponsored  by  the  Air  Force  Weapons  Laboratory  at  Kirtland  Air  Force  Base, 
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NM.  This  author  takes  no  credit  for  INVERT  but  did  check  its  function  to  ensure  it  was 
accurate  in  its  calculations. 

3. 4-2.3  MMAE  Subroutine.  The  subroutine  MMAE,  contained  in  entirety  in 
Appendix  A  in  Figure  A.l,  was  written  to  perform  nearly  all  of  the  MMAE-related  calcu¬ 
lations  in  MMSOFE.  The  argument  MOUT  which  is  used  in  calling  MMAE,  determines 
what  MMAE’s  function  will  be  as  follows: 

1.  Initialize  the  problem  at  the  start  of  the  simulation  (Calls  INITPR). 

2.  Initialize  the  runs.  INITPR  is  called,  several  parameter-related  variables  are  initial¬ 
ized,  and  the  probabilities  for  all  the  elemental  filters  are  initialized  to  one  divided 
by  the  number  of  elemental  filters. 

3.  Copy  the  vectors  saved  from  the  last  propagation  and  update  cycle  for  an  elemental 
filter  into  working  vectors  for  the  current  cycle.  The  variables  saved  are  as  follows: 

•  AFS  =  FILTER/SYSTEM  MAPPING  MATRIX 

•  XF  =  FILTER  STATE  VECTOR 

•  XS  =  SYSTEM  TRUTH  STATE  VECTOR 

•  EF  =  FILTER/SYSTEM  ERROR  VECTOR 

•  KF  =  FILTER  KALMAN  GAIN  VECTOR 

•  CF  =  FILTER  CONTINUOUS  FEEDBACK  CONTROL  MATRIX 

•  DF  =  FILTER  DISCRETE  FEEDBACK  CORRECTION  MATRIX 

•  FF  =  FILTER  FUNDAMENTAL  DYNAMICS  MATRIX 

•  QF  =  FILTER  PROCESS  NOISE  POWER  DENSITY  (STRENGTH)  MA¬ 
TRIX 

•  PC  =  COMPOSITE  COVARIANCL,  » ECTOR,  UPPER  TRIANGLE 

•  PF  =  FILTER  COVARIANCE  VECTOR,  UPPER  TRIANGLE 

•  PE  =  FILTER  ERROR  COVARIANCE  VECTOR,  UPPER  TRIANGLE 

•  YT  =  TRAJECTORY  VECTOR  AT  TIME  T 
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As  discussed  in  Section  3.4.4  with  regard  to  vectorization,  not  all  of  these  would  have 
strictly  required  saving  for  the  satellite  orbit  problem,  but  it  was  deemed  prudent  in 
order  to  make  future  adaptations  easier. 

4.  Save  the  vectors  shown  in  #3  directly  above  for  the  next  iteration. 

5.  Write  the  following  MMAE-related  values  to  output  file  designated  for  the  MMAE- 
related  data  (BLCNTOM): 

•  T  =  Simulation  time. 

•  XFMAT  =  State  estimates  for  all  elemental  filters 

•  XFM  =  Blended  (MMAE)  state  estimates. 

•  PFMM  =  Blended  state  covariance  matrix. 

•  PFMMD  =  Diagonals  of  the  blended-state  covariance  matrix. 

•  EM  =  Blended  state  estimation  error. 

•  PROB  =  Probabilities  of  the  elemental  filters. 

•  PAREST  =  Parameter  estimate. 

•  EA  =  Psirameter  estimation  error. 

•  EP  =  Diagonals  of  the  parameter  estimation  error. 

•  TRUPAR  =  True  value  of  the  parameter  from  subroutine  HRZSYS. 

6.  Save  xj,/,,,.  for  residual  calculations  and  save  Pjji,*,.  into  matrix  form  for  calculation 
of  the  covariance  of  the  blended  residuals. 

7.  Perform  the  following  MMAE  calculations  required  for  each  elemental  filter. 

•  Save  the  current  elemental  filter’s  state  estimates  and  upper  triangular  covari¬ 
ance  vector  values  into  the  storage  matrices,  XFMAT  and  PFMAT. 

•  Calculate  the  measurement  residual  vector. 

•  Calculate  the  covariance  matrix  (ALPHAM)  of  the  measurement  residual  vector 
and  find  the  inverse  of  ALPHAM. 
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•  Calculate  the  exponeatial  value  used  to  calculate  the  density. 

•  Calculate  the  coefficient  value  which  precedes  the  exponential  used  to  calculate 
the  density, 

8.  Perform  calculations  for  Bayesian  multiple  model  adaptive  estimation  only  when  the 
last  elemental  filter  has  been  complete  for  each  propagation/update  cycle. 

•  Calculate  the  denominator  of  the  conditional  probabilities. 

•  Evaluate  the  conditional  probabilities. 

•  Check  the  conditional  probabilities  against  the  minimum  probability  value  set 
in  the  —INPUT  files.  If  any  are  smaller: 

-  Set  them  equal  to  the  minimum  probability  value. 

-  Recalculate  the  denominator  of  the  new  adjusted  conditional  probabilities. 

-  Re-evaluate  the  conditional  probabilities. 

•  Calculate  the  blended  (MMAE)  state  estimates  using  the  probabilities  as  weight¬ 
ing  factor  on  the  individual  state  estimates  from  the  elemental  filters. 

•  Calculate  the  covariance  matrix  for  the  blended  state  estimates  (see  Equation 
(2.31)). 

•  Compute  the  blended  state  estimatation  error  from  the  true  state  values. 

•  Calculate  the  estimates  of  the  parameters. 

3.4. 2.4  OPFILE  Subroutine.  The  subroutine,  OPFILE,  was  written  to  open 
and  name  the  input  and  output  files  in  MMSOFE  automatically.  The  OPFILE  subroutine 
is  contained  in  its  entirety  in  Appendix  A,  Figure  A.2.  OPFILE  assigns  names  and  FOR¬ 
TRAN  unit  file  numbers  to  the  MMSOFE  input  and  output  files  as  required  by  the  number 
of  elemental  filters  being  utilized.  The  naming  convention  is  that  the  normal  ASCII  file¬ 
name  has  a  two  digit  numerical  prefix  corresponding  to  the  elemental  file  number.  An 
extra  few  files  beginning  with  “BL”  instead  of  the  elemental  filter  number  are  also  opened, 
with  “BL”  representing  “BLended”.  The  FORTRAN  unit  file  numbers  are  assigned  with 
an  even  hundred’s  value  representing  the  type  of  file  with  a  value  added  which  corresponds 
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to  the  elemental  filter  number,  or  one  more  than  the  total  number  of  elemental  filters  for 
the  “BL”  files.  For  example,  if  there  were  three  elemental  filters, 

•  The  MSOFE  CTOM  file  had  a  FORTRAN  file  unit  number  =  2 

•  The  MMSOFE  files  corresponding  to  CTOM  would  be  OlCNTOM,  02CNTOM, 
03CNTOM,  and  BLCNTOM. 

•  The  FORTRAN  unit  file  numbers  for  these  __CNTOM  files  would  be  201,  202,  203, 
and  204,  respectively. 

OPFILE  was  constructed  to  be  fiexible  to  the  extent  that,  if  the  computer  being  used  could 
handle  such  a  large  problem  quickly  enough,  OPFILE  would  allow  MMSOFE  to  process 
an  MMAE  simulation  with  up  to  98  elemental  filters  (i.e.  xx  =  98).  These  file  names  as 
corresponding  to  MSOFE  file  names  are  shown  in  Table  3.2  on  page  3-9,  and  the  vectorized 
variables  in  OPFILE  are  discussed  in  Section  3.4.4. 

3.4.3  Versions  of  MSOFE  Subroutines. 

3.4.3. 1  MAIN  Subroutine.  In  MMSOFE,  the  MAIN  program  is  no  longer  the 
top-level  executive  routine,  but  rather  that  distinction  has  been  relegated  to  SIMRUN.  In 
fact,  MSOFE’s  MAIN  subroutine  (see  Figure  3.1)  was  stripped  of  nearly  all  of  its  functions 
in  developing  MMSOFE  MAIN  (see  Figure  3.4).  In  comparing  the  two  versions  of  the 
MAIN  subroutine,  as  shown  in  these  Figures,  this  difference  is  apparent.  The  BLKDAT 
“call”  is  handled  as  an  “External”  call,  not  as  a  subroutine  call  in  both  versions.  The 
subroutine  calls  in  MSOFE  MAIN  to  subroutines  INIPR,  INITRN,  and  OUT  have  been 
moved  down  into  SIMRUN  and  the  new  MMSOFE  subroutine  (MMAE)  will  also  be  called 
from  SIMRUN,  The  only  true  subroutine  call  remaining  in  the  MMSOFE  MAIN  is  the  call 
to  SIMRUN,  with  a  loop  around  tue  “CALL  SIMRUN”  statement  to  facilitate  incrementing 
the  run  number  for  Monte  Carlo  simulations.  So,  in  comparing  Figure  3.1  (MSOFE  MAIN) 
to  Figure  3.4  (MMSOFE  MAIN),  the  INITPR,  INITRN,  and  OUT  subroutines  and  the 
one  new  subroutine  MMAE  are  the  only  changes  to  the  boxed  list  of  subroutines  that  are 
caUed  by  SIMRUN. 
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No 


Figure  3.4.  MMSOFE  MAIN  Program 
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3. 4- 3. 2  SIMRUN  Subroutine.  Beside  being  used  as  the  MMSOFE  top-level 
executive  subroutine,  SIMRUN  is  also  required  to  handle  the  coordination  of  which  elemen¬ 
tal  filter  is  currently  being  processed.  The  program  flow  for  the  new  MMSOFE  SIMRUN 
subroutine  is  shown  in  Figures  3.5  and  3.6.  This  is  accomplished  primarily  via  calls  to 
the  new  MMAE  subroutine,  monitoring  the  elemental  filter  number,  and  using  the  flags, 
BGNRUN  and  BGNPROB.  BGNPROB  stays  “true”  until  all  elemental  filters  are  initial¬ 
ized  at  the  start  of  a  simulation.  BGNRUN  remains  “true”  while  initializing  all  elemental 
filters  at  the  start  of  each  Monte  Carlo  run.  Comparing  Figure  3.2  to  Figures  3.5  and  3.6, 
other  changes  from  MSOFE  to  MMSOFE  include; 

1.  MMSOFE  includes  time  manipulations  to  ensure  a  common  time  progression  for  all 
elemental  filters. 

2.  In  MMSOFE,  the  ENDRUN  flag  is  set  to  “true”  only  when  the  last  elemental  filter 
has  finished  the  last  run. 

3.  The  calls  to  subroutine  OUT  are  controlled  both  during  a  run  and  at  the  end, 
OUT(8),  of  the  simulation  for  each  elemental  filter  in  MMSOFE.  Subroutine  0UT(8) 
was  called  from  MAIN  in  MSOFE. 

4.  MMSOFE’s  SIMRUN  decrements  the  run  counter  (IRUN)  between  elemental  filters 
during  run  initializations.  In  both  MSOFE  and  MMSOFE,  IRUN  is  incremented  in 
subroutine  INITRN  and  therefore,  it  must  be  decremented  between  i)asses  through 
INITRN  for  each  elemental  filter. 

Therefore,  it  should  be  obvious  by  comparing  the  SIMRUN  subroutine  from  MSOFE  versus 
SIMRUN  from  MMSOFE  that  the  basic  structure  of  MSOFE  remains  intact  in  MMSOFE. 

3. 4- 3.3  MEASUP  Subroutine.  The  MEASUP  modifications  can  be  seen  by 
comparing  MSOFE  MEASUP  depicted  in  Figure  3.3  with  MMSOFE  MEASUP  as  shown  in 
Figure  3.7.  MEASUP  in  MMSOFE  makes  three  calls  to  subroutine  MMAE.  The  additions 
to  MSOFE  MEASUP  are  the  three  calls  to  the  MMAE  subroutine,  some  initializations  at 
the  beginning  of  MEASUP,  and  the  vectorizations  discussed  in  Section  3.4.4. 
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Figure  3.5.  MMSOFE  SIMRUN  Subroutine  (Part  A) 
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Figure  3.6.  MMSOFE  SIMRUN  Subroutine  (Part  B) 
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Figure  3.7.  MMSOFE  MEASUP  Subroutine 
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Once  again,  the  changes  to  MSOFE  are  fairly  minimal,  preserving  the  MSOFE  struc¬ 
ture  and  using  it  for  calculating  each  independent  elemental  filter’s  values  in  turn.  Because 
the  changes  were  minimal,  no  diagram  of  subroutines  HRZSYS  nor  HRZFIL  was  made, 
but  those  changes  were  important  from  an  MMAE  point  of  view.  Besides  inserting  the 
failure  signals  in  HRZSYS  and  HRZFIL,  the  respective  truth  and  filter  measurements  or 
measurement  noise  variances  were  loaded  into  matrix  form  for  use  in  MMAE  calculations 
in  subroutine  MMAE.  Additionally,  in  HRZFIL,  the  filter  derivative  vector,  HFILF,  was 
also  saved  into  a  matrix  for  MMAE  calculation  use.  Because  the  calculations  were  done  for 
each  elemental  filter  independently  and  in  turn  during  a  propagation/measurement  update 
cycle,  these  matrices  were  re-used  for  the  elemental  filters  instead  of  having  to  be  saved 
for  each  elemental  filter  until  MMAE  calculations  could  be  accomplished. 

3.4. 3.4  MM ATH  Subroutines.  The  MMATH  subroutine  is  merely  a  very  use¬ 
ful  linear  algebra  utility  subroutine  which  simplifies  calculations  which  utilize  vectors  and 
matrices.  No  diagram  of  the  program  flow  for  the  MMATH  subroutine  was  included  herein 
either  because,  in  large  part,  MMATH  was  duplicated  from  CSOFE.f.  There  were  some 
additions  made  to  MMATH  for  use  in  MMSOFE.  In  comparison  to  the  original  MMATH, 
the  version  of  MMATH  used  in  MMSOFE  has  the  following  added  options: 

1.  New  option  to  find  the  determinant  of  a  square  matrix  of  any  size. 

2.  New  option  to  transform  a  vector  containing  the  upper  triangular  terms  of  a  matrix 
into  a  symmetric  matrix.  While  this  operation  and  the  next  operation  are  performed 
in  MSOFE  via  specific  FORTRAN  code,  they  are  not  performed  in  the  MMATH 
subroutine  in  MSOFE. 

3.  New  option  to  create  a  vector  containing  the  upper  triangular  terms  of  a  symmetric 
matrix. 

4.  Option  to  copy  one  column  of  an  array  to  one  column  of  another  array,  which  can 
be  used  for  building  matrices  from  vectors  or  for  pulling  one  column  vector  out  of  a 
matrix. 
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In  order  to  clarify  the  MMAE  code  (included  in  Figure  A.l  in  Appendix  A),  which 
contains  a  number  of  calls  to  subroutine  MMATH,  a  brief  description  of  the  MMATH  call 
statement  is  given. 


CAU  IOUTH(  R.RR.CR.  A.RA.CA.  B.RB.CB  ) 


The  ‘‘R,RR,CR,”  “A,RA,CA,’’  and  ‘‘B,RB,CB,”  refer  to  the  resultant  array  (R)  and 
the  two  input  arrays  (A  and  B)  with  the  number  of  rows  (RR,  RA,  and  RB)  and  number 
of  columns  (CR,  CA,  and  CB)  following  the  name  of  each  array.  The  matrix  functions 
possible  with  this  version  of  MMATH  are: 
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>»> 

MATRIX  EQUIVALENCE 

R  * 
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»♦» 

MATRIX  ADDITION 
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R(, 

)  *  AO 
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*■  A(,) 
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COL  Rf.I)  <-  COL  A(,J) 

R(,I)  =  A(,J) 

'D' 

MATRIX  DETERMINATE 
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<A> 

'1' 

"VECTOR"  INNER  PRODUCT 
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A~  *  B 

'0' 

"VECTOR"  OUTER  PRODUCT 
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'C' 

VECTOR  CROSS  PRODUCT 

R  = 

A  X  B 

'N' 

VECTOR  NORM  (MAGNITUDE) 

R  = 

<A> 

'U' 

UNIT  VECTOR 

R  * 

A  /  <A> 

Several  of  these  need  some  added  explanation. 


•  “MATRIX  FILL”  =  A  scalar  value  passed  in  a  1  X  1  matrix  is  loaded  into  all  elements 
of  matrix  R. 

•  “COL  R(,I)  <-  COL  A(,J)”  =  The  “J”  column  of  matrix  A  is  copied  into  the  “I” 
column  of  matrix  R.  This  can  be  used  for  stripping  off  one  column  vector  from  a 
matrix  composed  of  column  vectors. 
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3.4-4  Vectorization  and  Matricization.  Once  OPFILE  was  completed,  the  “vector- 
ization”  of  the  appropriate  variables  in  MSOFE.f  and  USOFE.f  was  a  major  endeavor 
required  before  multiple-filter  simulations  of  any  kind  could  be  done.  The  number  of 
variables  which  were  vectorized  was  kept  to  a  minimum,  although  some  variables  were 
vectorized  because  it  would  be  necessary  for  other  applications,  even  though  it  might  not 
have  been  necessary  for  this  problem.  The  variables  which  needed  to  be  vectorized  were 
declared  in  either  the  24  COMMON  BLOCKS  or  in  the  subroutine-contained  local  decla¬ 
ration  statements.  These  variables,  as  existing  in  MSOFE.f,  are  discussed  in  the  MSOFE 
User’s  Manual  (24)  and  are  well  described  in  MSOFE.f  and  USOFE.f  programs  themselves. 
Therefore,  they  will  only  be  identified  as  having  been  vectorized,  but  their  functions  won’t 
be  described,  except  if  vitally  relevant. 

Figure  3.8  shows  the  common  blocks  which  were  created  for  MMSOFE  that  did  not 
exist  in  MSOFE.  In  Figure  3.8,  the  variable  name  is  given  first  with  its  dimensions  followed 
by  it  data  type  (IN  =  INTEGER,  RL  =  REAL,  DP  =  DOUBLE  PRECISION,  and  etc.). 
Next  is  given  a  function  descriptor  such  that 

•  “U”  =  Used  in  the  subroutine. 

•  “S”  =  Set  in  the  subroutine. 

•  “A”  =  Argument  passed  in  the  CALL  statement. 

COMMON  BLOCK  MMVAR  contains  almost  all  of  the  MMAE-related  variables.  While 
many  of  the  variables  in  MMVAR  could  have  been  locally  declared  as  local  variables  in 
subroutine  MMAE,  they  were  declared  as  COMMON  variables  in  MMVAR  to  ease  future 
adaptations  utilizing  the  MMAE-related  variables.  With  these  variables  in  COMMON,  a 
subroutine  for  another  application  could  be  added  to  USOFE.f  and  implemented  merely 
by  adding  a  call  statement  in  MMSOFE,  since  all  of  the  MMAE-related  variables  are 
available  through  the  MMVAR  COMMON.  Further,  no  effort  was  made  economize  on  the 
number  of  variables,  but  instead  a  new  variable  was  utilized  for  each  new  MMAE-related 
value  calculated.  This  makes  it  easy  to  follow  the  calculations  through  the  code,  and  it 
also  makes  it  easier  to  write  intermediate  variables  out  to  a  files  during  development.  The 
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MMVAR  variables  are  described  in  Figure  3.9,  where  they  are  arranged  in  logical  order, 
instead  of  alphabetically. 

The  COMMON  BLOCK,  DBUG,  is  vital  to  MMSOFE  multiple  model  operation. 
It  contains  TTELF,  which  is  used  as  the  maximum  number  of  files  throughout  MSOFE, 
ELF  (elemental  filter)  which  is  the  integer  corresponding  to  the  elemental  filter  currently 
being  utilized  and  the  xum_  vectors  containing  the  FORTRAN  unit  numbers  referred  to 
in  Table  3.2.  The  number  suffixes  on  the  xum_  vectors  correspond  to  the  hundreds  digit 
of  the  MMSOFE  unit  number  or  the  MSOFE  unit  number  for  the  corresponding  sets  of 
files.  If  a  “BL"  file  is  accessed,  either  ELFT  (ELF  temporary)  or  TTELF  are  used,  not 
ELF. 

The  vectorized  variables  and  some  new  associated  variables  are  specified  in  Figures 
3.10  and  3.11  in  accordance  with  where  they  are  declared.  In  these  figures  and  in  MMSOFE, 
the  variable,  TELF,  equals  the  total  number  of  elemental  filters,  while  the  variable  TTELF 
and  the  parameter  DUMI  both  equal  TELF  +  1.  The  reason  for  having  two  is  because 
DUMI  is  a  parameter  and  is  used  to  dimensionalize  vectors  in  subroutine  OPFILE;  DUMI 
cannot  be  modified  in  the  program,  whereas  TTELF  is  the  true  number  of  each  type  of  file 
to  open  in  subroutine  OPFILE.  Therefore,  TTELF  =  TELF  +  1  (the  number  of  elemental 
filters  plus  one  for  the  “BL”  files)  unless  TELF  =  1,  in  which  case,  TTELF  =  1  because 
no  “BL”  files  are  opened.  The  MMSOFE  filenaming  convention  was  discussed  in  Section 
3.4.  It  is  apparent  from  an  inspection  of  Figure  3.11  that  INIT  (local  initialization  flag) 
was  the  only  local  variable  vectorized  in  USOFE.f,  albeit  in  five  USOFE.f  subroutines,  and 
INIT  was  vectorized  in  most  of  the  subroutines  in  MMSOFE.f.  Figure  3.11  also  lists  all 
of  the  local  variables  in  MMSOFE.f  which  were  vectorized.  Several  special-use  vectors 
were  also  declared  as  local  variables  for  the  simple  purpose  of  using  them  as  arguments  in 
subroutine  call  statements.  The  only  reason  for  the  concern  is  that,  if  a  call  statement  to 
a  subroutine  in  MSOFE  had  a  vector  for  an  argument  which  was  matricized  in  MMSOFE, 
the  subroutine  might  not  be  passing  the  row  or  column  desired.  On  some  computers  or 
with  some  FORTRAN  compilers,  matrices  are  stored  in  memory  according  to  columns, 
while  on  other  computers,  the  storage  could  be  by  rows.  This  method  of  using  the  special- 
use  vectors  obviated  this  argument  ambiguity  and  eliminated  any  possibility  of  error  in 
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passing  between  subroutines  one  row  or  one  column  of  a  matrix  (in  MMSOFE)  which  had 
formerly  been  a  vector  (in  MSOFE).  This  practice  also  ensures  the  compatability  of  the 
code  with  other  computer  systems  and  FORTRAN  compilers. 


COMDECK  MMVAR 

COMMON/  MMVAR  /  XFMAT,HFILFM,ZFM,ZSM,RSM,RFM,PFMAT,XFM, 
PFMNS,XFMNS,HXFMNS,ZRESM,ZRESMT,HPF, 
HFILFMT,  HPHTF,ALPHAM.ALPHAMI,ALPHAMIR, 
RTAR,DET,PROB.  DPROB.DXF.PFCOL.PFM.PFMPLS, 
MINPROB,PFMM,DXFT,ZM,EM,DXFMAT,PFMMV,YY, 
PFMMVT,PFMMD,PAREST,PARMN,TRUPAR,EA,EP, 
EXPROB.EXCOEF.RTARV.PARCNT,  MINFLAG 
REAL  XFMAT(NF,TELF),  HFILFM(NZ,NF),  ZFM(NZ,1), 

ZSM(NZ).  RSM(NZ,NZ),  RFM(NZ.NZ). 

PFMAT{NTF,TELF),XFM(NF).  PFMNS(NF,NF), 

XFMNS(NF,1).  HXFMNS(NZ,1),  ZRESM(NZ,1). 

ZRESMT(1,NZ),  HPF(NZ,NF),  HFILFMT(NF,NZ), 

HPHTF(NZ,NZ),  ALPHAM(NZ,NZ),  ALPHAMI(NZ,NZ), 
ALPHAMIR(NZ,1).  RTAR,  DET(NZ,NZ), 

PROB(TELF,l),  DPROB,  DXF(NF), 

PFCOL(NTF),  PFM(NF,NF),  PFMPLS(NF,NF), 

MINPROB,  PFMM(NF,NF).  DXFT(1,NF), 

ZM(NZ),  EM(NF),  DXFMAT(NF,TELF), 

PFMMV(NTF),  YY(NIN),  PFMMVT(NTF), 

PFMMD(NF),  PAREST(DUMI.PN),PARMN(PN), 

TRUPAR(PN),  EA(PN),  EP(PN), 

EXPROB(TELF,l),  EXCOEF(TELF),  RTARV(TELF) 

INTEGER  PARCNT,  MINFLAG 
COMDECK  DBUG 

COMMON/  DBUG  /  TTELF.ELF 

XUM02,  XUM03,  XUM04,XUM05,  XUM06,  XUM07, 

XUM08,  XUM09,  XUMlO,  XUM12,  XUM15 
INTEGER  TTELF.ELF 

XUM02(DUMI),  XUM03.  XUM04(DUMI),XUM05(DUMI), 
XUM06(DUMI),  XUM07(DUMI),XUM08(DUMI), 

XUM09,  XUMlO,  XUM12,  XUMl5(DUMI) 

COMDECK  DBUG2 

INTEGER  DBG 


Figure  3.8.  NEW  COMMON  BLOCK  Variables 
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ARGUMENTS 

MOUT  IN  U  A 

COMMON  VARIABLES  ACCESSED 
/MMVAR/ 

XFMAT(NF,TELF)  RL()  US 
HFILFM(NZ,NF)  RL()  US 
HFILFMT(NF,NZ)  RL()  US 
ZFM(NZ,1)  RL()  US 

ZSM(NZ)  RL()  US 

RFM(NZ.N4j  RL()  US 
RSM(NZ,NZ)  RL()  US 
PFMAT(NTF,TELF)RL()  US 
XFM(NF)  RL()  US 

PFMNS(NF.NF)  RL()  US 
XFMNS(NF,1)  RL()  US 
HXFMNS(NZ,1)  RL()  US 

ZRESM(NZ,1)  RL()  US 
ZRESMT(l.NZ)  RL()  US 
HPF(NZ,NF)  RL()  US 

HPHTF(NZ,NZ)  RL()  US 
ALPHAM(NZ,NZ)  RL()  US 
ALPHAMI(NZ,NZ)  RL()  US 
ALPHAMIR(NZ.l)  RL()  US 
RTAR  RL()  US 

DET(NZ,NZ)  RL()  US 
PROB(TELF,l)  RL()  US 
DPROB  RL()  US 

DXF(NF)  RLO  US 

DXFT(1,NF)  RLO  US 
DXFMAT(NF,TELF)RLO  US 
PFM(NF,NF)  RL()  US 
PFCOL(NTF)  RLO  US 
MINPROB  RLO  U 

PFMM(NF,NF)  RLO  US 
PFMMVT(NTF)  RLO  US 

PFMMV(NTF)  RLO  US 

EM(NF)  RLO  US 

PAREST(DUMI,PN)RLO  US 
TRUPAR(PN)  RL()  US 
EA(PN)  RLO  US 

EP(PN)  RLO  US 

EXPROB(TELF,l)  RLO  US 
EXCOEF(TELF)  RL()  US 
RTARV(TELF)  RLO  S 
MINFLAG  IN  US 

PN  IN  U 

NZ  IN  U 


INDICATES  ACTION  TO  TAKE 


MATRIX  OF  XF  -  ALL  ELEMENTAL  FILTERS 

HFILF  IN  A  2D  MATRIX 

TRANSPOSE  OF  HFILM 

VECTOR  OF  FILTER  MEASUREMENTS 

VECTOR  OF  SYSTEM  MEASUREMENTS 

VECTOR  OF  FILTER  MEASUREMENT  NOISE  VARIANCES 

VECTOR  OF  SYSTEM  MEASUREMENT  NOISE  VARIANCES 

MATRIX  OF  PF  UPPER  DIAG  -  ALL  ELEMENTAL  FILTERS 

MMAE  BAYESIAN  BLENDED  STATE  ESTIMATES 

PF  MINUS 

XF  MINUS 

VECTOR  OF  FILTER  MEASUREMENTS  FOR  EL 
RESIDUAL  CALCULATION  HXFMNS  =  HFILFM  *  XFMNS 
EL  FILTER  RESIDUAL 
TRANSPOSE  OF  ZRESM 
HFILFM  *  PFMNS 

HFILFM  *  PFMNS  *  HFILFM(TRANSPOSE) 

COVARIANCE  OF  MEASUREMENTS,  R+H*P*H(TRANSPOSE) 
ALPHAM(INVERSE) 

ALPHAM(INVERSE)  *  RESIDUAL 

.5  ♦  RESIDUAL(TRANSPOSE)*ALPHA(INVERSE)*RESIDUAL 
DETERMINANT  OF  THE  ALPHAM 

VECTOR  OF  CONDITIONAL  PROBABILITY  FOR  EL  FILTERS 
DENOMINATOR  OF  CONDITIONAL  PROBABILITY 
EL  FILTERS  -  BLENDED  STATE  ESTIMATES  (XF-XM) 
DXF(TRANSPOSE) 

MATRIX  OF  DXF’S  FOR  ALL  EL  FILTERS 

COVARIANCE  MATRIX  W/  DXF  (P=DXF*DXF(TRANSPOSE) 

COLUMN  OF  UPPER  TRIANGLE  OF  PFM 

MINIMUM  PROBABILITY  BOUND  (SET  IN  INPUT  FILES) 

COVARIANCE  OF  BLENDED  ESTIMATES 

COVARIANCE  OF  EACH  ELEMENTAL  FILTER 

PFMMVT(I)=PROB(ELFT,l)*(PFMAT(I,ELFT)+PFCOL(I)) 

BLENDED  COVARIANCE  PFMMV=PROB(ELF)*PFMMV(ELF) 

PFMMV(I)=PFMMV(I)+PFMMVT(I) 

BLENDED  ESTIMATION  ERROR  (EM  =  XFM  -  AFS  *  XS) 

BAYESIAN  ESTIMATE  OF  THE  PARAMETERS 

TRUE  VALUE  OF  THE  PARAMETER  FROM  HRZSYS 

ERROR  IN  THE  PARAMETER  ESTIMATE 

DIAGONAL  OF  PARAMETER  CALCULATED  COVARIANCE 

EL  FILTER  DENSITY  (EXCOEF  *  EXP(RTAR)) 

LEADING  COEFFICIENT  OF  DENSITY  FOR  EL  FILTER 
VECTOR  OF  RTAR  FOR  ALL  EL  FILTERS,  (OUTPUT  ONLY) 
FLAG  SET  IF  MINPROB  IS  USED 
NUMBER  OF  MEASUREMENTS 
NUMBER  OF  MEASUREMENTS 


Figure  3.9.  NEW  COMMON  BLOCK  Variables  Descriptions 
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COMDECK  CORBIT 

COMMON  /CORBIT/  ZBOUND,  ZFBIAS,  RFNOIS, 

NWRF,  NWRS,  NWTF,  NWTS,  RNERR, 

RANGLF,  RRANGF,  RANGLS,  RRANGS 
REAL  TMPULS(2)MU,RHO,ZBOUND(TELF),ZFBIAS(TELF) 

RFNOIS(TELF),NWRF(TELF),  NWRS(TELF),NWTF(TELF) 
NWTS(TELF),RANGLF(TELF),RRANGF(TELF), RANGLS, RRANGS 
INTEGER  RNERR(TELF) 

COMDECK  SIZES 

INTEGER  TELF,XUN15,DUMI,PN,NZ 

PARAMETER  (PN  =  2),  (XUNI5  =  101),  (TELF  =  5),(DUMI  =  TELF  +  1), 

(TNF  =  TELF*NF),  (TNS  =  TELF*NS),(NZ  =  2) 

COMDECK  XUNITS 

INTEGER  XUN02,XUN03,XUN04,XUN05,XUN06,XUN07,XUN08,XUN09,XUN10, 
XUN12,XUN15,XUM02(DUMI),XUM03,XUM04(DUMI),XUM05(DUMI), 
XUM06(DUMI),XUM07(DUMI),XUM08(DUMI),XUM09,XUM10,XUM12 
PARAMETER  (XUN02=201,  XUN03=3,  XUN04=401,  XUN05=501,  XUN06=601) 
(XUN07=701,  XUN08=801,  XUN09=9,  XUN10=10,XUN12=12) 

COMDECK  EVCNTL 

COMMON/  EVCNTL  /  DTFDBK,  DTNOYS 
REAL  DTFDBK(TELF),  DTNOYS(TELF) 

COMDECK  ICBLK 

COMMON/  ICBLK  /  DV 
REAL  DV(TELF,2) 

COMDECK  INCOND 

COMMON/  INCOND  SEDXFI,  SEDXSI 
INTEGER  SEDXFI(TELF),  SEDXSI(TELF) 

COMDECK  KOUNTS 

COMMON/  KOUNTS  /  NPP 
INTEGER  NPP(TELF) 

COMDECK  MESDO 

COMMON/  MESDOA  /  ALPHA,  ZF,  ZRES,  ZS,  ZTOL 
REAL  ALPHA(TELF),ZF(TELF),ZRES(TELF),ZS(TELF),ZTOL(TELF) 
COMDECK  MODFIL 

COMMON/  MODFIL  /  HF,  HFILF,  HFILS 

REAL  HF(TELF,NF),HFILF(TELF,NF),  HFILS(TELF,NS) 

COMDECK  MODSYS 

COMMON/  MODSYS  /  HS,  HSYSF,  HSYSS 

REAL  HS(TELF,NS),  HSYSF(TELF,NF),  HSYSS(TELF,NS) 

COMDECK  MONOFF 

COMMON/  MONOFF  /  LCOV 
LOGICAL  LCOV(TELF) 


Figure  3.10.  COMMON  BLOCK  variables 
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USOFE.f  Local  VuiaUes 


DECK  DFDBAK 

INTEGER  INIT(TELF) 

DECK  FCQFIL 

INTEGER  INIT(TELF) 

DECK  FCQSYS 

INTEGER  INIT(TELF) 

DECK  FSMAP 

INTEGER  INIT(TELF) 

DECK  XFDBAK 

LOGICAL  INIT(TELF) 

MMSOFE.f  Local  Variablea 

DECK  DERIV 

INTEGER  INIT(TELF) 

REAL  XSFDTI(TELF,NS+NF),DXSFDTI(  NF+NS  ) 
DECK  EPRINT 

LOGICAL  INIT(TELF) 

DECK  GETP 

INTEGER  TELF 
DECK  KUTMER 

INTEGER  INIT(TELF) 

REAL  HDYN(TELF) 

LOGICAL  ATSTEP(TELF) ,  ATTOUT(TELF) 

DECK  07PL0T 

INTEGER  INIT(TELF),ELFT 
INTEGER  NXPTS(TELF) 

DECK  OUT 

INTEGER  MACEPT(TELF) ,  MREJCT(TELF) 
DECK  RDGAIN 

INTEGER  INIT(TELF) 

DECK  TNEXTE 

LOGICAL  INIT(TELF) 

DECK  TNEXTM 

INTEGER  INIT(TELF) 

DECK  VALDAT 

LOGICAL  DATAOK(TELF) 

DECK  OPFILE 

INTEGER  ELFT.DUMI.TTELF.I 
CHARACTER*?  FNAME(DUMI) 

CHARACTER*2  FNUMA(DUMI) 

CHARACTER*!  FNUMIOO(DUMI).  FNUMlO(DUMI) 
CHARACTER*!  FNUM!(DUM1) 

INTEGER  FNUM(DUMl) 

DECK  MMAE 

REAL  XSO(TELF,NS),  XFO(TELF,NF) 

REAL  PCO(TELF,NTC),  PFO(TELF,NTF) 

REAL  PEO(TELF,NTF),  YTO(TELF,NYT) 

REAL  EFO(TELF,NF),  AFSO{TELF,NF,NS) 

REAL  CFO(TELF,NF,NF),  DFO(TELF,NF,NF) 
REAL  FFO(TELF,NF,NF),  QFO(TELF.NF,NF) 
REAL  KFO(TELF,NF) 


Figure  3.11.  MMSOFE.f  and  USOFE.f  Local  Variables 
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The  vectorized  local  variables  in  the  two  new  subroutines,  the  OPFILE  and  MMAE 
subroutines,  also  require  some  explanation.  In  OPFILE,  FNAME  is  the  vector  containing 
the  entire  file-names,  including  the  elemental  filter  number  or  “BL”  as  derived  by  OPFILE. 
Similarly,  the  vector,  FNUM,  contains  the  fORTRAN  unit  number  associated  with  a  fiie. 
ELFIOO,  ELFIO,  and  ELFl  each  correspond  to  the  integer  value  of  the  hundreds,  tens,  and 
ones  digit  of  the  unit  number,  while  the  same  representation  in  character  form  is  given  in 
variables  FNUMlOO,  FNUMIO,  and  FNUMl.  These  are  used  to  derive  the  prefix  on  the 
file  names  which  correspond  to  the  elemental  filter  number  or  “BL”  for  blended  (see  Table 
3.2). 

Because  MMSOFE  propagates  and  updates  each  elemental  filter  in  turn  before  con¬ 
tinuing  to  a  new  propagation/update  cycle  for  any  of  the  filters,  vectors  are  required  to 
save  the  current  elemental-filter-specific  values  until  the  propagations  and  updates  have 
been  accomplished  for  all  of  the  elemental  filters.  The  saved  values  for  an  elemental  filter 
are  then  reloaded  when  it  is  that  filter’s  turn  to  propagate  again.  This  permits  the  inde¬ 
pendent  but  parallel  processing  of  the  xx  elemental  filters.  These  local  vectors  are  declared 
in  the  MMAE  subroutine  as  shown  in  Figure  3.11.  The  MMSOFE-pertinent  COMMON 
BLOCK  variables  are  shown  in  Figure  3.10.  Two  variables  whose  values  are  supplied  via 
the_INPUT  files  were  shown  in  Figure  3.10  just  to  emphasize  that  they  weren’t  vectorized. 
These  vectors  are; 

1.  TMPULSE  =  times  to  apply  thrusts  for  Hohmann  transfer. 

2.  RHO  =  radius  ratio,  apogee/perigee. 

They  weren’t  vectorized  because  these  are  part  of  the  true  (system)  orbit  determinations 
and  the  true  orbit  is  the  same  for  all  elemental  filters.  They  could  be  easily  vectorized 
if  required  for  some  specific  application.  The  other  problem-specific  (Group  5)  variables 
from  the  .JNPUT  files  which  were  vectorized  are: 

1.  NWRF  =  filter  range- rate  process  noise  strength  (note  “R”  for  “range”). 

2.  NWRS  =  system  range-rate  process  noise  strength. 

3.  NWTF  =  filter  angle-rate  process  noise  strength  (note  “T”  for  “theta”). 

4.  NWTS  =  system  angle-rate  process  noise  strength. 
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5.  RRANGF  =  filter  range  measurement  noise  variance. 

6.  RRANGS  =  system  range  measurement  noise  variance. 

7.  RANGLF  =  filter  angle  measurement  noise  variance. 

8.  RANGLS  =  system  angle  measurement  noise  variance. 

9.  ZBOUND  =  bound  on  acceptable  measurement  residual  size,  in  multiple  of  filter- 
computed  residual  standard  deviation,  beyond  which  a  single  measurement  will  be 
edited  out  (no  update  performed). 

10.  MINPROB  =  lower  bound  on  elemental  filter  probabilities. 

11.  RNERR  =  Number  to  specify  the  failure  type  to  be  simulated. 

12.  ZFBIAS  =  filter  measurement  bias. 

13.  RFNOIS  =  scale  factor  on  filter  measurement  noise. 

With  some  of  these  variables,  vectorization  may  not  have  been  strictly  necessary,  but  it 
was  implemented  in  order  to  simplify  modifications  later  on.  Essentially,  if  there  was  a 
question  of  whether  to  vectorize  or  not  to  vectorize  at  this  time,  it  was  chosen  to  “err”  and 
vectorize,  rather  than  be  ultra-conservative  and  make  future  adaptations  more  difficult 
than  necessary. 

All  of  the  “XUN0_”  parameters  in  sizes  and  XUNITS  provide  the  starting  number  for 
OPFILE  to  use  in  generating  the  input  and  output  files.  In  the  cases  of  XUN03,  XUN09, 
XUNIO,  and  XUN12,  these  are  files  correspond  to  one  single  MMSOFE  file  (see  Table  3.2). 
All  of  the  other  vectorized  variables  given  in  Figure  3.10  required  vectorization  in  order  to 
keep  the  data  and  calculations  separate  for  the  elemental  filters. 

3.5  MMSOFE  Validation 

The  MMSOFE  validation  corresponds  to  those  software  validation  tests  described  in 
Section  1.6.4  as  follows: 

1.  The  failures  chosen  for  MMSOFE  validation  (see  Section  1.6.3  and  Table  1.1)  were 
programmed  into  subroutine  HRZSYS  to  facilitate  using  the  variable  RNERR  in  the 
—INPUT  files  to  designate  the  failure  type. 
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2.  The  failure  types  shown  in  Table  1.1  were  run  with  MSOFE  to  verify  the  tuning 
parameters  and  hlter  performance. 

3.  The  MMSOFE  software  was  partially  validated  using  multiple,  but  identical  elemen¬ 
tal  filters,  without  the  implementation  of  MMAE  probability  and  blending  being 
computed.  The  MMSOFE-generated  elemental  filter  data  were  compared  to  ensure 
all  elemental  filters  were  performing  identically  and  then  compared  to  data  obtained 
from  corresponding  MSOFE  runs.  These  data  sets  were  compared  “bit-for-bit.” 
This  comparison  was  done  using  the  Matrix*  { 10)  software  to  difference  a  large  num¬ 
ber  of  pertinent  variables  and  then  plot  those  differences.  An  exact  plot  of  zero  for 
the  differences  indicated  success.  These  numerous  plots  of  straight  lines  at  zero  are 
not  included  in  this  thesis  for  obvious  reasons.  Suffice  it  merely  to  list  the  variables 
which  were  compared  and  checked.  These  variables  were: 

•  XF  =  Filter  state  estimates 

•  XS  =  System  (truth)  state  estimates 

•  EF  =  Filter/System  state  estimation  error 

•  ZF  =  Filter  (predicted)  measurement 

•  ZS  =  System  (actual)  measurement 

•  PF  =  Vector  of  the  upper  triangular  of  the  filter  covariance  matrix 

•  PC  =  Vector  of  the  upper  triangular  of  the  composite  covariance  matrix 

•  PE  =  Vector  of  the  upper  triangular  of  the  filter  error  covariance  matrix 

•  ZRES  =  Measurement  residual  (actual  -  estimated) 

•  ALPHA  =  Estimated  variance  of  the  measurement  residual  which  is  calculated 
in  MSOFE  as  a  scalar  iterative  update  rather  than  being  calculated  from  a 
ZRES  vector.  This  means  there  cannot  be  any  “off  diagonal”  terms  in  this 
ALPHA. 

4.  The  same  values  were  again  compared  to  validate  the  MMSOFE  software  by  using 
multiple,  but  distinct  elemental  filters  and  comparing  the  elemental  filter  results  to 
data  obtained  from  matching  MSOFE  runs.  Again  the  zero  plots  are  not  included. 
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5.  Again  the  values  were  compared  after  the  Bayesian-blending  MMAE  calculations 
were  included.  These  comparisons  were  accomplished  both  with  and  without  a 
Hohmann  transfer  in  the  profile. 

6.  The  blended  state  estimates  were  compared  to  the  state  estimates  obtained  from  the 
elemental  filters  and  “sanity-checked”  against  plots  of  the  residuals  and  probabilities 
generated  by  the  elemental  filters.  Additionally,  comparisons  of  the  calculations  were 
made  to  verify  the  accuracy  of  the  Bayesian-blending  code  in  the  subroutine  MMAE. 
Several  passes  through  the  multiple  model  MMAE  calculations  were  compared  with 
identical  calculations  done  using  the  Matrix^  (10)  software.  This  comparison  was 
done  for  calculations  both  before  and  after  the  time  at  wh'ch  the  failures  were  induced 
and  repeated  for  several  different  failure  types.  No  comparison  plots  were  made  of 
these  comparisons,  but  the  numerical  values  were  verified  visually  on  the  computer 
screen  to  an  accuracy  of  nine  significant  digits. 

3.6  Summary 

The  goals  of  MMSOFE  development  have  been  to  develop  a  tool  for  simulating 
multiple  model  adaptive  estimation  simulations  which  would  be  robust  and  easy  to  use  for 
anyone  familiar  with  MSOFE.  Additionally  and  to  test  the  new  software,  MMSOFE  should 
be  able  to  indicate  failures  via  observing  the  Bayesian  probabilities  and  also  providing  an 
estimate  of  the  uncertain  parameters.  These  goals  have  been  met,  and  the  MMSOFE  tool 
developed  should  be  as  transportable  to  other  computer  systems  as  MSOFE. 
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IV.  Analysis  of  Results 


4.1  Overview 

This  chapter  discusses  the  results  of  the  MMSOFE  multiple  model  simulations  for 
the  five  failure  types  discussed  in  Section  1.6.3  and  shown  in  Table  1.1,  plus  two  simu¬ 
lations  which  were  added  to  those  in  Table  1.1.  Table  4.1  shows  all  seven  failure  types, 
including  Cases  #6  and  #7,  which  were  added  to  help  clarify  unforeseen  MMAE-related 
phenomenon  described  in  Section  4.3.2.  Cases  #6  and  #7  and  the  rationale  for  adding 
them  are  described  in  Section  4.4.2.  The  plots  which  present  the  data  from  the  simula¬ 
tions  for  each  of  the  seven  failure  types  are  located  in  Appendix  B,  grouped  into  sections  of 
plots  corresponding  to  each  fsdlure  type.  Preceding  each  section  of  plots  is  also  the  related 
portion  of  Table  4.1.  The  explanations  and  variable  dehnations  located  with  Table  1.1  in 
Section  1.6.3  also  apply  to  Table  4.1.  The  plots  are  scaled  uniformly  according  to  what 
data  is  presented,  with  the  plot  scale  chosen  so  as  to  include  all  of  the  data  even  in  the 
worst  data  set.  In  other  words,  that  all  of  the  probability  plots  have  a  scale  from  0  to  1  and 
the  residual  plots  and  the  state  estimation  error  plots  all  have  a  scale  between  +  .5  and  - 
.5.  For  Cases  #1,  #2,  and  #3,  the  parameter  error  plots  are  scaled  from  -  .01  to  -|-  .01  and 
the  parameter  estimate  plot  range  is  from  -.002  to  .01.  For  Cases  #4,  #5,  #6,  and  #7,  the 
parameter  error  plots  are  scaled  from  -  .2  to  -f  .2  and  the  parameter  estimate  plot  range 
is  from  -.03  to  .15.  By  standardizing  the  plot  scales  as  much  as  possible  in  this  manner, 
it  should  be  easier  to  make  good  comparisons  between  data  presented  on  different  plots. 
The  plots  corresponding  to  each  failure  type  are  organized  as  indicated  in  the  introduction 
of  Appendix  B  and  also  as  follows: 

1.  Figure  B.l:  One  plot  of  the  probability  associated  with  each  elemental  filter. 

2.  Figure  B.2  and  Figure  B.3:  Figure  B.2  range  and  Figure  B.3  for  angle,  each  figure 
containing  five  separate  plots.  Each  separate  plot  is  related  to  an  elemental  filter 
and  graphs  the  data  for: 

•  The  mean  residual. 

•  Mean  residual  ±  Monte  Carlo  calculated  standard  deviations. 
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Table  4.1.  Orbit  Problem  Failure  Cases 


Error 
Case  # 

Failure  Type 

In  Truth  Model 

Failure 

Induced 

Elemental  Filter 

1 

Step  in 

Measurement  Noise 

Rs  —  Rso  +  fi-Rso 

T  >  25 

R/o(l  + 

2 

Ramp  in 

Measurement  Noise 

25  <  T  <  35 

T  >  35 

R/o(l  +  BIASfij ) 

3 

Step/Return  in 
Measurement  Noise 

Rs  —  Rso 

Rs  =  Rso  +  6Rso 

T  <  25,  T  >  35 
25  <  T  <  35 

R/o(l  +  BIAShj) 

4 

Step  in 
Measurement 

T  >  10 

Zj-\-  BI AS^f 

5 

Ramp  in 
Measurement 

10  <  T  <  20 

T>  20 

Zf  -|-  BIASzf 

6 

Step  in 
Measurement 

T  >  10 

Zj  BI ASzf 

7 

Ramp  in 
Measurement 

10  <  T  <  25 

T  >  25 

Zj  -h  BIASzj 

•  Zero  mean  db  filter-predicted  standard  deviations. 

3.  Figures  B.4  to  B.7:  Four  figures,  one  corresponding  to  each  of  the  four  states  (range, 
range-rate,  theta,  and  theta-rate).  Each  figure  consists  of  six  separate  plots.  The 
plots  correspond  to  the  blended  and  five  elemental  filters  and  depict  data  for: 

•  The  mean  state  estimation  error. 

•  Mean  estimation  error  ±  Monte  Carlo  ccilculated  standard  deviations. 

•  Zero  mean  ±  filter-predicted  standard  deviations. 

4.  Figure  B.8:  One  figure  with  two  plots  relating  to  the  MMAE/Bayesian  parameter 
estimation: 

•  One  plot  of  the  mean  estimation  error,  mean  estimation  error  ±  Monte  Carlo 
standard  deviation,  and  zero  ±  filter-predicted  standard  deviations. 

•  The  second  plot  of  the  true  parameter  and  mean  estimated  parameter  values. 

Plots  occur  sequentially  thereafter  in  such  groups  of  eight. 
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4-2  Elemental  Filter  Description 

The  five  elemental  filters  (#1,  #2,  #3,  #4,  #5)  were  designed  with  Filter  #3  always 
being  the  baseline  filter  -  no  biases  applied  to  the  measurement  nor  to  the  measurement 
noise  variance.  For  all  of  the  simulations: 

•  Filter  #1  has  a  measurement  bias  which  matched  the  final  value  of  the  truth  bias 
magnitude  (except  for  Case  #3  of  Table  4.1). 

•  Filter  #2  has  a  bias  equal  to  |  of  Filter  #l’s. 

•  For  Cases  #4  to  #7  which  have  measurement  bias/ramp-type  failures,  Filters  ^4 
and  #5  had  measurement  biases  opposite  in  sign  but  equal  in  magnitude  to  those  of 
Filters  #2  and  #1,  respectively. 

•  For  Cases  #1  to  #3  with  the  failure  induced  via  a  bias/ramp  added  to  the  mea¬ 
surement  noise  variance,  negative  values  for  the  measurement  noise  variance  would 
not  be  correct,  so  the  Filter  variance  bias  halved  the  baseline  noise  variance 
magnitude  and  Filter  #5  halved  Filter  #4’s  variance  bias  magnitude.  Referring  to 
Table  4.1,  this  halving  was  accomplished  by  setting  BIASr^  =  -.5  in  Case  #4  and 
BIASrj  —  -.75  in  Case  #5,  which  corresponds  to  the  labels  on  Figures  B.l  through 
B.23. 

The  particular  discretization  of  the  parameter  space  utilized  in  this  thesis  made  it  easier 
to  observe  probability  results  but  it  is  not  necessarily  the  best  discretization  for  state  or 
parameter  estimates.  Refinement  of  the  parameter  discretization  was  not  a  consideration 
in  this  thesis  effort. 

4.3  Probabilities 

4-3.1  Expected  Probabilities.  Because  Filter  model  4^3  matches  the  truth  model  in 
the  absence  of  failures,  one  might  expect  its  probability  to  be  highest  when  there  is  no 
failure.  Additionally,  Filter  44^1  matches  the  true  situation  when  the  induced  failure  has 
reached  its  maximum  (full-failure)  and  Filter  #2  is  midway  between  Filters  #3  and  ^1 
and  therefore  Filter  ^2  matches  the  situation  midway  between  no-failure  and  full-failure 
(important  particularly  in  ramp-type  failure  modes).  Heuristically,  it  might  be  expected 


that  the  probabilities  in  Case  #2  (ramping  of  the  measurement  noise  variance)  should 
shift  from  Filter  #3,  pass  to  Filter  #2  as  the  failure  magnitude  grows,  and  then  shift 
over  to  Filter  #1  when  the  failure  is  fully  culminated.  In  Cases  #1  and  #3,  it  might  be 
expected  that  the  probabilities  would  switch  directly  from  Filter  #3  to  Filter  #1  with 
only  a  small  increase  (if  any)  in  Filter  #2’s  probability.  These  expectations  were  proven 
to  be  essentially  valid  for  Cases  #1,  and  #3  in  which  the  failure  was  induced  in  the 
measurement  noise  variance  (see  Figures  B.l,  B.9,  B.17),  but  not  for  the  other  cases  in 
which  the  failure  was  induced  directly  via  the  measurement  (see  Figures  B.25  and  B.33). 
Filter  #2’s  probability  did  increase  after  failure  induction  for  Cases  #1  to  #3,  but  Filter 
^2  never  became  dominant. 

4-3.2  Unexpected  Probability  Phenomena.  The  heuristics  broke  down  somewhat  for 
the  cases  in  which  the  failure  was  induced  as  a  “weak”  ramp  directly  on  the  measurements 
(Cases  #5,  See  Table  4.1),  growing  to  a  bias  equal  in  magnitude  (.1  distance  units)  to  the 
measurement  failure  bias  (Cases  ^^14).  They  also  broke  down  for  Case  #4  in  which  the  bias 
was  introduced  suddenly  as  a  “small”  step  (Appendix  B,  Section  B.4).  Both  the  “small” 
step  and  the  “weaJc”  ramp  cases  (Appendix  B,  Section  B.5)  culminated  with  a  bias  of  .1 
distance  units  in  the  measurement.  In  both  of  these  cases,  the  probabilities  did  not  shift 
all  the  way  over  to  Filter  #1  but  stayed  “hung-up”  on  the  intermediate  Filter  #2.  At 
least  they  shifted  away  from  the  baseline  filter  (#3),  but  sticking  on  Filter  #2  was  not  the 
expected  result. 

This  phenomenon  can  be  understood,  though,  by  examining  the  residuals  which  are 
the  sources  of  the  information  for  the  hypothesis  conditional  probability  calculations  within 
the  MMAE  and  by  keeping  in  mind  that  a  probability  determined  at  one  measurement 
update  time  carries  over  into  the  probability  calculation  performed  at  the  next  update 
time  by  direct  weighting  of  the  density  function  (Equation  (2.27)).  Therefore,  once  Filter 
#2  dominated  the  probability,  its  residuals  would  have  to  become  very  “bad”  and  another 
filter’s  residuals  would  have  to  become  “good”  to  shift  probability  from  Filter  #2  to  that 
other  filter.  Recognition  of  this  unexpected  phenomenon  is  made  more  important  because 
of  the  relationship  of  these  cases  to  spoofing.  Ail  four  cases,  #4,  #5,  #6,  and  #7  (Section 
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4.4.2),  represent  the  simulation  of  spoofing  which  was  discussed  in  Chapter  1  and  by 
Vasquez  (36)  with  regard  to  a  Global  Positioning  System.  This  phenomenon  indicates  how 
difficult  detection  and  isolation  of  gradually  growing  spoofing-type  errors  really  are,  even 
when  using  MMAE  techniques. 

4.4  Analysis  of  Simulation  Data 

4-4-1  Failure  Induced  by  Measurement  Noise  Variance  Bias  or  Ramp.  Cases  #1, 
#2,  and  #3  (see  Figures  B.l  to  B.24)  all  represent  failures  akin  to  signal  jamming  in 
the  Global  Positioning  System,  as  discussed  in  Chapter  1  and  by  Vasquez  (36).  This  is 
a  simpler  failure  type  to  detect  than  the  spoofing  type  failures  of  Cases  #4  through  #7. 
Consequently,  there  were  no  great  surprises,  as  those  which  will  be  described  in  Section 
4.4.2.  Case  #1  represents  the  total  loss  of  GPS  signal  or  signal  jamming  after  T  =  25. 
Case  #2  represents  a  gradual  fading  out  of  the  signal  or  gradual  application  of  a  jamming 
signal  from  T  =  25  to  T  =  35  and  then  maintaining  of  the  level  at  the  end  of  the  ramp  for 
the  remainder  of  the  simulation.  Case  #3  represents  the  same  conditions  as  Case  #1  but 
with  the  jamming  or  signal  loss  ending  at  T  —  35. 

Initially,  during  the  benign  (no-failure)  part  of  all  three  of  these  simultations,  the 
probabilities  for  Filters  #2,  #4  and  baseline  all  vied  for  the  highest  probability.  It  seems 
heuristically  obvious  that  these  three  filters  should  have  the  highest  probabilities  because 
these  three  filters  most  closely  match  the  truth  in  the  no-failure  situation.  By  about 
T  =  12,  the  baseline  elemental  filter  had  become  dominant  with  the  highest  probability 
and  remained  dominant.  The  fact  that  the  baseline  filter  ultimately  remained  dominant 
in  the  absence  of  any  f<ulure  was  confirmed  by  one  simulation  with  Tjinai  =  100  (the  plots 
of  which  are  not  included  in  this  thesis). 

In  Cases  #1  and  #3  which  had  the  failure  induced  via  a  measurement  noise  variance 
bias  beginning  at  T  =  25,  the  probability  began  shifting  from  the  baseline  filter  to  Filter  #1 
immediately  after  inception  of  the  failure.  With  Case  #3,  the  probability  began  shifting 
back  from  Filter  #1  to  the  baseline  filter  when  the  failure  was  discontinued  at  T  =  35. 
The  probability  shifted  much  more  slowly  from  the  baseline  filter  to  Filter  #1  in  Case  #2 
because  the  ramped  failure  caused  a  much  slower  change  in  the  residuals  and  therefore  a 
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slower  shifting  of  the  probability.  Figure  B.IO  is  key  to  understanding  why  the  probability 
mass  leaves  Filter  #3  and  is  attracted  to  Filter  #1  vs.  sequentially  to  Filter  #2  and  then 
to  Filter  #1.  By  the  time  Filter  #3’s  residuals  are  poor  enough  to  cause  ps  to  decrease, 
Filter  #l’s  residuals  are  better  (relative  to  the  filter-computed  Ai)  than  those  of  Filter 
#2.  Plot  B.16  shows  good  identification  of  Filter  #2,  though  it  is  very  similar  to  Plot  B.8 
in  which  the  parameter  estimate  would  be  anticipated  to  converge  more  quickly  to  Oi- 

It  is  also  interesting  to  note  the  following  similarities  in  the  probabilities  calculated 
for  Cases  #1,  #2,  and  #3. 

1.  The  probabilities  for  Cases  #1,  #2,  and  #3  were  identical  (as  they  should  be)  for 
all  five  filters  until  T  =  25  when  the  failures  were  induced. 

2.  They  were  identical  until  T  =  35,  for  Cases  #1  and  #3,  at  which  time  the  failure  in 
Case  #3  was  turned  off  (see  Figures  B.l  and  B.17). 

3.  Likewise,  the  probabilities  for  the  five  filters  were  almost  the  same  for  Cases  #1  and 
#2  by  the  end  of  the  simulation  at  T  =  50  because  the  failure  offset  had  been  the 
same  for  these  two  Cases  since  T  =  35. 

4.  Presumably,  the  probabilities  for  Cases  #1  and  #2  would  have  become  equal  some¬ 
time  after  T  =  50. 

4-4-^  Failure  Induced  by  Measurement  Bias  or  Ramp.  Referring  to  Table  4.1  and 
to  the  figures  in  Appendix  A,  which  depict  the  residuals  from  the  measurement  ramp  and 
bias  failure  cases,  the  following  observations  are  of  interest; 

1.  The  angle  residual  plots  (Figures  B.27,  B.35,  B.43  and  B.51)  provided  no  clear  in¬ 
sights. 

2.  As  anticipated,  the  baseline  filter  had  the  best  range  residuals  and  the  highest  prob¬ 
abilities  before  the  failure  onset  in  all  four  of  these  cases  (#4,  #5,  #6,  and  #7). 

3.  For  Cases  #4  and  #5  (.1  bias  and  ramp  cases): 

•  Filter  #2  attracted  and  kept  the  highest  probability  from  the  baseline  filter  right 
after  the  failure  was  induced  and  while  Filter  #1,  which  intuitively  should  have 
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acquired  the  highest  probability,  responded  slightly  in  probability  accumulation 
but  then  lost  out  to  Filter  #2. 

•  Examination  of  the  corresponding  range  residuals  (Figures  B.26  and  B.34)  re¬ 
veals  that  the  Filter  #2  range  residual  (mean  residual  values  over  the  15  Monte 
Carlo  runs)  was  within  one  filter-predicted  standard  deviation  of  zero  through¬ 
out  the  entire  simulation.  The  mean  range  residual  for  Filter  #1  was  outside  one 
filter-predicted  standard  deviation  of  zero  before  failure  induction  even  though 
it  was  essentially  zero  after  the  full  f<ulure  realization.  It  was  hypothesized  that 
the  probability  could  be  forced  toward  the  “correct”  filter  either  by  allowing 
a  longer  simulation  time  over  which  Filter  #l’s  residuals  would  have  better 
characteristics  than  those  in  Filter  #2,  or  by  letting  the  truth  model  maximum 
offset  be  larger  so  as  to  make  Filter  #2's  residuals  more  distinctly  poor  relative 
to  Filter  ^I’s  residuals.  The  hypothesis  of  allowing  the  longer  simulation  time 
was  tried  without  the  desired  result.  This  then  led  to  Cases  #6  and  #7. 

4.  For  Cases  #6  and  #7  (.15  bias  and  ramp  cases): 

•  While  Filter  #2  began  to  attract  the  probability  from  the  baseline  filter  just  as 
in  Cases  #4  and  #5,  Filter  #1  soon  took  over  and  had  the  highest  probability 
for  the  rest  of  the  simulation,  just  as  originally  hypothesized. 

•  Examining  the  corresponding  range  residuals  clarifies  the  situation.  While  the 
residuals  are  identical  to  those  of  Cases  #4  and  #5  before  induction  of  the 
fciilure,  the  larger  magnitude  of  the  final  bias  value  caused  the  range  residuals 
of  Filter  #1  to  be  worse  than  they  were  in  Cases  #4  and  #5,  but  caused  those 
for  Filter  #2  to  become  much  worse.  In  fact,  the  range  residuals  for  Filter  #2 
are  slightly  more  than  one  standard  deviation  from  zero  and  those  for  Filter  #1 
somewhat  less.  Thus,  the  probability  can  switch  from  the  baseline  filter  first  to 
Filter  #2  and  finally  over  to  Filter  #1,  as  originally  anticipated. 

These  observations  support  the  postulate  presented  in  Section  4.3.2.  The  range  residuals 
from  Cases  #4  and  #5  were  better  for  Filter  #2,  and  the  Filter  #2  probability  higher 
than  those  for  Filter  #1,  immediately  after  the  failure  was  introduced.  Therefore,  Filter 
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#2  continued  to  attain  higher  and  higher  probabilities  according  to  Equation  (2.26).  By 
the  time  Filter  #rs  residuals  became  better  than  Filter  #2’s,  Filter  #2  had  good  residuals 
and  higher  probabilities,  so  Filter  #1  couldn’t  attract  the  probability  over  the  time  of  the 
simulation.  However,  with  Cases  #6  and  #7,  Filter  #2’s  range  residuals  became  so  much 
worse  than  those  of  Filter  #2,  that  Filter  #2  gave  up  the  probability  to  Filter  #1,  which 
had  better  range  residuals. 

The  observations  with  regard  to  Cases  #4  and  #5  were  the  impetus  for  the  last  two 
test  cases  which  hadn’t  originally  been  planned  but  were  retro-fitted  into  Table  1.1  as  the 
test  Cases  #6  and  #7  in  Table  4.1.  Case  #6  has  a  truth  measurement  failure  bias  of  .15 
distance  units  in  contrast  to  .1  distance  units  of  Case  #4  (see  Appendix  B,  Section  B.6  and 
Figures  B.41  -  B.48).  Case  #7  is  the  case  in  which  the  failure  was  induced  as  a  “strong” 
ramp  on  the  measurement  (see  Appendix  B,  Section  B.7  and  Figures  B.49  -  B.56),  The 
failure  was  induced  identically  to  the  prior  “weak”  ramp  test  case  (Case  #5),  with  the 
ramps  of  equal  slope  but  with  ramp  endurance  of  1  j  times  as  long  for  Case  #5.  Therefore, 
for  Case  #7,  the  failure  bias  on  the  measurement  was  .15  instead  of  just  .1  distance  units 
after  the  slope  ended  at  T  =  25. 

4-4 -3  State  Estimation.  The  MMAE-blended  state  estimates  for  all  four  states 
appear  quite  accurate  from  the  state  estimation  error  plots  shown  in  Appendix  B.  In  each 
of  these  cases,  the  blended  estimation  error  mean  value,  standard  deviation  of  that  error, 
and  the  filter-predicted  standard  deviation  were  smaller  than  they  were  for  any  of  the  five 
elemental  filters.  While  not  accomplished  in  this  research,  a  good  performance  baseline 
could  be  created  by  running  a  single  filter  simulation  in  which  the  filter  is  artificially 
informed  of  the  exact  nature  and  time  of  failure. 

The  accuracy  of  the  state  estimates  for  Cases  #4  -  #7  was  excellent  oefore  the  failures 
were  induced,  with  mean  estimation  errors  of  nearly  zero  and  small  standard  deviations  of 
the  estimation  error.  After  the  failures  were  induced,  the  blended  state  estimation  error 
was  as  small  as  the  state  estimation  error  corresponding  to  the  elemental  filter  with  the 
highest  probability.  Thus,  for  Cases  #4  and  #5,  the  errors  were  as  small  as  Filter  #2’s 
(highest  probabiUty)  even  though  Filter  #1  had  smaller  state  estimation  errors,  essentially 
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zero  mean.  After  introduction  of  the  failure,  the  standard  deviations  of  the  estimation 
errors  followed  a  pattern  similar  to  that  of  the  state  estimation  error  mean  values  which 
could  be  expected  because  they  were  probability  dependent  also. 

4-4’4  Parameter  Estimation.  Before  the  failures  were  induced  in  the  simulations, 
the  parameter  estimates  were  nearly  exact  for  all  of  the  cases  simulated.  However,  after 
the  failure  was  induced,  because  Filter  #1  in  each  of  the  first  five  error  cases  modelled 
the  final  magnitude  of  the  induced  failure  exactly,  the  parameter  estimates  obtained  in 
all  of  these  five  cases  (Table  4.1)  must  be  smaller  than  the  value  hypothesized  in  Filter 
^1.  Additionally,  the  parameter  estimate  could  only  be  that  large  if  the  probability  for 
Filter  #1  were  equal  to  1,  which  is  not  possible  since  the  other  probabilities  were  lower 
bounded  at  .01.  This  upper  bounding  of  the  parameter  estimate  is  due  to  the  fa.ct  that 
the  MMAE-based  parameter  estimate  is  the  sum  of  the  probability-weighted  parameters 
hypothesized  in  each  elemental  filter.  If  the  estimated  parameter  were  equal  to  or  greater 
than  the  true  value,  the  probabilities  would  also  have  to  sum  to  something  greater  than 
one.  That  is  not  a  possible  solution,  and  thus  the  parameter  estimates  must  be  biased  in 
this  application. 

The  parameter  estimation  errors  for  Cases  #1  -  #3  were  generally  smaller  in  the 
steady  state  after  the  failure  was  fully  realized  than  for  Cases  #4  -  #7.  For  Cases  #4  - 
#7,  the  steady  state  parameter  estimation  error  mean  values  were  all  nearly  equal  (-.05). 
For  Cases  #6  and  #7  in  which  the  hypothesized  parameter  value  in  Filter  #1  was  .1  and 
the  true  parameter  value  was  .15,  worse  parameter  estimates  might  have  been  anticipated. 
However,  except  for  the  transient  period  during  and  slightly  after  the  failure  was  applied, 
the  parameter  estimation  error  was  nearly  the  same  as  for  Cases  #4  and  ^5. 

Failure  Detection  and  Isolation  (FDl)  with  MMAE  (see  Sections  1.5.4  and  2.3)  is 
accomplished  merely  by  monitoring  the  probabilities  for  each  of  the  elemental  filters  or, 
equivalently,  by  monitoring  the  parameter  estimate.  For  the  test  cases  (see  Table  4.1) 
used  in  this  thesis,  the  elemental  filter  which  modelled  the  correct  failure  was  detected  for 
all  cases  except  Cases  ^4  and  #5.  In  these  two  cases.  Filter  #2  was  identified  instead 
of  correctly  identifying  Filter  #1.  This  difficulty  could  possibly  have  been  overcome  by 
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additional  discretization  studies  in  which  more  elemental  filters,  some  with  hypothesized 
parameters  exceeding  the  possible  true  parameter  values.  That  might  have  resulted  in 
correct  FDI  and  improved  the  state  and  parameter  estimates  as  well.  Another  possibility 
to  overcome  the  difficulty  is  to  restart  the  probabilities  periodically  at  the  lower  bound 
value  for  all  elemental  filters  except  the  baseline  filter  (perhaps  in  parallel  with  the  original 
MMAE  probabilities).  In  this  scenario,  as  probability  “dumps  out  of”  Filter  #3  when 
its  residuals  are  large,  the  probability  mass  ought  to  be  absorbed  by  the  filter  with  the 
best-looking  residuals. 

4.5  Summary. 

Chapter  IV  has  described  the  results  of  the  MMSOFE  orbit  problem  simulations. 
Analysis  of  the  elemental  filter  performance,  the  blended  state  estimates,  the  probabilities 
and  parameter  estimatation,  and  the  FDI  results  were  discussed. 


4-10 


V.  Conclusions  and  Recommendations 


This  chapter  provides  a  brief  discussion  of  conclusions  reached  with  regard  to  the 
MMSOFE  software  and  the  simulation  results.  Recommendations  for  future  research  are 
also  presented. 

5. 1  Conclusions 

5.1.1  MMSOFE.  The  primary  goal  of  this  thesis  (see  Section  1.2)  was  to  develop 
a  simulation  and  an<ilysis  software  tool  to  aid  controls  designers  and  system  integrators 
in  developing  and  analyzing  MMAE  applications.  As  described  in  Chapter  III,  the  input 
and  output  file  organization  used  by  MSOFE  (24)  was  preserved.  As  discussed  in  detail  in 
Chapter  IV,  MMSOFE  performed  the  Kalman  filter  calculations  for  each  elemental  filter  of 
the  MMAE,  which  was  verified  by  bit-for-bit  data  comparisons  with  data  produced  using 
the  same  filter  models  in  MSOFE.  MMSOFE  also  calculates  the  probabilities  associated 
with  the  correctness  of  each  elemental  filter,  calculates  the  MMAE  blended  state  estimates 
and  MMAE-computed  state  estimation  error  covariance,  and  estimated  the  parameter  val¬ 
ues.  MMSOFE  will  process  from  one  to  ninety-eight  elemental  filters  merely  by  modifying 
the  file  SIZES  and  by  generating  an  input  file  for  each  elemental  filter. 

Inherent  to  the  primary  goal  were  several  other  goals: 

•  The  MSOFE  structure  was  modified  as  little  as  possible  in  developing  MMSOFE. 

•  The  MMSOFE  software  tool  should  be  user-friendly  for  anyone  already  familiar  with 
MSOFE 

•  The  method  of  implementing  MMAE  via  MMSOFE  could  make  it  easy  to  adapt 
MMSOFE  for  other  multiple  model  applications  as  well.  The  method  of  using  the 
subroutine  SIMRUN  as  the  MMSOFE  top-level-executive  routine  and  to  handle  the 
coordination  of  which  elemental  filter  is  being  processed  makes  the  MMSOFE  adap¬ 
tation  for  other  applications  quite  simple.  Thus  MMSOFE  could  be  easily  expanded 
to  accommodate  applications  such  as  MMAE-based  control,  multiple  model  adaptive 
control,  or  distributed  Kalman  filtering. 
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5.1.2  Simulations.  The  detailed  discussion  of  the  results  from  the  satellite  orbit 
estimation  problem  was  presented  in  Chapter  IV  and  will  not  be  repeated.  Suffice  it  to 
say  that: 

•  The  performance  of  MMSOFE’s  calculations  for  each  elemental  filter  was  verified  as 
accurate  based  on  comparison  to  MSOFE. 

•  The  MMAE-based  calculation  of  probabilities,  state  estimates,  parameter  estimates, 
and  FDI  for  the  cases  simulated  were  entirely  reasonable  and  were  deemed  accurate. 

5.2  Recommendations 

A  brief  list  of  recommendations  for  future  research  is  presented,  with  many  of  the 
details  concerning  these  items  presented  earlier. 

•  Adapt  MMSOFE  for  MMAE-based  control. 

•  Adapt  MMSOFE  for  multiple  model  adaptive  control. 

•  Adapt  MMSOFE  for  distributed  Kalman  filtering  and  compare  the  results  with  the 
results  from  other  DKF  software  tools,  such  as  Distributed  Kalman  Filter  Simulation 
(DKFSIM)  (3). 

•  Use  MMSOFE  to  apply  MMAE  to  the  integrated  navigation  reference  system  in¬ 
vestigated  by  Vasquez  in  his  thesis  (36).  The  navigation  reference  system  integrates 
an  inertial  navigation  system,  a  range/range-rate  transponder  system,  and  a  Global 
Positioning  System. 

•  Use  MMSOFE  for  other  applications,  such  as  those  listed  in  Section  1.5.3. 

•  Since  MMSOFE  estimates  the  parameters  using  probabilty  weighting,  MMSOFE 
could  be  used  as  an  aid  to  tune  a  filter  optimally.  This  could  be  accomplished  by 
distributing  the  discretized  parameter  space  across  the  elemental  filters  as  with  any 
MMAE  simulation,  letting  the  uncertain  parameters  be  the  filter  tuning  parameters. 
The  estimate  of  the  parameter  could  then  be  used  for  subsequent  MMSOFE  or 
MSOFE  simulations. 

•  The  previous  recommendation  could  be  changed  slightly  and  MMSOFE  modified 
to  take  partial  derivatives  of  variables  of  interest  with  respect  to  the  discretized 


parameter  values.  These  variables  of  interest  could  be  state  estimates,  covariances, 
or  estimation  errors  (or  their  statistics).  These  partial  derivatives  could  then  be 
utilized  by  a  gradient  algorithm  to  help  find  optimal  tuning  parameters  or,  for  that 
matter,  to  help  find  optimal  values  for  uncertain  parameters  located  in  any  part  of 
the  filter  model. 
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Appendix  A.  New  MMSOFE  Subroutines 


This  appendix  contains  the  MMSOFE  FORTRAN  code  for  the  two  completely  new 
subroutines,  MMAE  and  OPFILE. 


A.l  MMAE  Subroutine 


*DECK  miAE 

SUBROUTIIE  mUE  (  NCHJT  ) 

C 

Cli.1  I  I— — — -  I .  —■.■I  —■■■I  ■  . . III,,.,., 

C  AVIOIICS  DIRECTORATE  NSOFE  l.S 

C  VRIGHT  LABORATORT  Jt3IE  1991 

C  VRIGHT-PATTERSOI  AFB  OHIO  COPYRIGHT 


C  HOTICE:  THIS  ROUTIHE  IS  SUBJECT  TO  THE  DISTRIBOTIOI  LINITATIOI. 
C  EXPORT  VARIIHG.  AID  IXSTROCTIOH  HOTICE  AT  TOP  OF  FILE. 


C  OUT  —  OUTPUT  COHTROLLER 
C 

USAGE  "  CALL  MMAE  (  MOOT  ) 

ABSTRACT  —  MMAE  COHTROLS  CQITROLS  THE  CALCUUTIOIS  REUTED  TO 
MULTIPLE  MODEL  CALCULATIOIS .  THE  TABLE  BELOV  GIVES 
C  MORE  DETAIL  ABOUT  HHAT  HAFPEIS  II  THE  IIDIVIDUAL  OUTPUT 

C  ROUniES  FOR  EACH  VALUE  OF  MOOT. 

C 

C 

C  LIMITATIOIS  ~  lOIE  KIOVI 
C 

C  REMARKS  ~ 

C  TABLE  OF  OUTPUT  ACTIOHSC**)  FOR  EACH  POSSIBLE  EVERT 

C  - 

C  MOOT  EVEIT 

C  IIDEX  DESCRIPTER 

C  - 

C  1  IIITIALIZE  PROBLEM 

C  2  IIITIALIZE  RUI 

C  3  COPT  VECTORS  TO  HORKIIG  VECTORS  PROM  LAST  ITERATIOI 

C  4  SAVE  V0RXII6  VECTORS  FOR  lEXT  ITERATIOI 

C  5  WRITE  MMAE  DATA  TO  OUTPUT  FILE  BLCITOM 

C  6  SAVE  XF  FOR  MMAE  STATE  ESTIMATES  A  SAVE  PF  IHTO  MATRIX  FORM 

C  FOR  CALCULATIOI  OF  REAL  ALFA 

C  7  PERFORM  MMAE  CALCUUTIOI  REQUIRED  FOR  EACH  ELEMEI7AL  FILTER 

C  8  PERFORM  CALCULATIOIS  FOR  MULTIPLE  MODEL  ADAPTIVE  ESTINATIOI 

C  OILY  VHEI  THE  LAST  ELEMEITAL  FILTER  HAS  BEEI  COMPLETE  FOR 

C  THIS  PROPAGATIOI/UPDATE  CYCLE 

C  - 

C  HEADER  DATE  —  93/10/29 


1 
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n  o 


C  AROONEITS 

C  HOOT  II  0  A  IIDICATES  ACTIOI  TO  TAKE 
C 

C  CCHOIOI  VARIABLES  ACCESSED 


c 

/mVAR/ 

c 

XFNAT(IF,TELF) 

RLO 

08 

MATRIX  OF  XF  -  ALL  ELEMEITAL  FILTERS 

c 

HFILFNdZ.IF) 

RLO 

US 

HFILF  II  A  2D  MATRIX 

c 

HFILFRT(IF.IZ) 

RLO 

US 

TRAISPOSE  OF  HFILM 

c 

ZFNdZ.l) 

RLO 

US 

VECTOR  OF  FILTER  MEASUREMEITS 

c 

ZSNdZ) 

RLO 

US 

VECTOR  OF  STSTEM  MEASUREMEITS 

c 

RFNdZ.IZ) 

RLO 

US 

VECTOR  OF  FILTER  MEASOREMEIT  lOISE  VARIAICES 

c 

RSHdZ.IZ) 

RLO 

OS 

VECTOR  OF  STSTEM  MEASOREMEIT  lOISE  VARIAICES 

c 

PFIUTdTF.TELF) 

RLO 

OS 

MATRIX  OF  PF  UPPER  DIAG  -  ALL  ELEMEITAL  FILTERS 

c 

XFNdF) 

RLO 

US 

MMAE  BATE8IAI  BLEIDED  STATE  ESTIMATES 

c 

PFnSdF.IF) 

RLO 

OS 

PF  MHOS 

c 

XFNISdF.l) 

RLO 

US 

XF  MUDS 

c 

HXFMISdZ.l) 

RLO 

US 

VECTOR  OF  FILTER  MEASUREMEITS  FOR  EL 

c 

RESIDUAL  CALCOUTIQI  HXFMIS  -  HFILFM  •  XFMIS 

c 

ZRESMdZ.l) 

RLO 

US 

EL  FILTER  RESIDUAL 

c 

ZRESNT(I.IZ) 

RLO 

US 

TRAISPOSE  OF  ZRESM 

c 

HPFdZ.IF) 

RLO 

US 

HFILFM  *  PFMIS 

c 

HPRTFdZ.IZ) 

RLO 

US 

HFILFM  *  PFMIS  *  HFILFM(TRaSPOSE) 

c 

ALPHANdZ.IZ) 

RLO 

OS 

COVARIAICE  OF  MEASUREMEITS.  R+H*P*H (TRAISPOSE) 

c 

ALPHAKIdZ.IZ) 

RLO 

US 

ALPHAMdlVERSE) 

c 

ALPHAKERdZ,!) 

RLO 

US 

ALPHAMdIVERSE)  *  RESIDUAL 

c 

RTAR 

RLO 

US 

.6  •  RESIDUAL(TRAISPOSE)*ALPHA(IIVERSE)»RESIOUAL 

c 

DETdZ.IZ) 

RLO 

US 

I»TERMIIAIT  OF  THE  ALPHAM 

c 

PROB(TELF.l) 

RLO 

US 

VECTOR  OF  COIDITIOIAL  PROBABILITT  FOR  EL  FILTERS 

c 

DPROB 

RLO 

OS 

DEIOMIIATOR  OF  COIDITIOIAL  PROBABILITT 

c 

OXFdF) 

RLO 

US 

EL  FILTERS  -  BLEIDED  STATE  ESTIMATES  (XF-XM) 

c 

DIFTCI.IF) 

RLO 

US 

DXF(TRAISPOSE) 

c 

DXFHATdF.TELF) 

RLO 

US 

MATRIX  OF  DXF’S  FOR  AU  EL  FILTERS 

c 

PFHdF.IF) 

RLO 

OS 

COVARIAICE  MATRIX  V/  DXF  (P-DXF*DXF  (TRAISPOSE) 

c 

PFCOLdTF) 

RLO 

US 

COLUm  OF  UPPER  TRIAIGLE  OF  PFM 

c 

MllPROB 

RLO 

u 

MIIIMOM  PROBABILITT  BODID  (SET  II  IIFDT  FILES) 

c 

PFNMdF.IF) 

RLO 

US 

COVARIAICE  OF  BLEIDED  ESTIMATES 

c 

PFMHVTdTF) 

RLO 

US 

COVARIAICE  OF  EACH  ELEMEITAL  FILTER 

c 

PFMMVT(I)-PROB(ELFT.  1)*(PFMAT(I  .ELFT)+PFCOL(I)) 

c 

PFmVdTF) 

RLO 

OS 

BIEIDED  COVARIAICE  PFMMV-PROB(ELF)*PFHMV(ELF) 

c 

PFMMV(I)-PFMRV(I)-^PFMRVT(I) 

c 

ENdF) 

RLO 

US 

BUEIDED  ESTIMATIOI  ERROR  (EM  -  XFM  -  AFS  *  XS) 

c 

BAREST (DOMI .PI) 

RLO 

US 

BATESIAI  ESTIMATE  OF  THE  PARAMETERS 

c 

TRDPAR(PI) 

RLO 

US 

c 

EACPI) 

RLO 

OS 

c 

EP(PI) 

RLO 

US 

c 

EXPROBCTELF.l) 

RLO 

US 

EL  FILTER  DEISITT  (EXCOEF  *  EXP (RTAR)) 

c 

EXCOEF(TELF) 

RLO 

OS 

LEADIIG  COEFFICIEIT  OF  DEISITT  FOR  EL  FILTER 

c 

RTARV(TELF) 

RLO 

s 

VECTOR  OF  RTAR  FOR  ALL  EL  FILTERS.  (OUTPUT  OILT) 

c 

NIIFU6 

11 

OS 

FLAG  SET  IF  MllPROB  IS  USED 

c 

/RDimi/ 

c 

IROIS  II  U 

lUNBER  OF  ROIS  II  SIMDUTIOI  STDDT 

c 

/OSRDTI/ 

c 

PC  RLO  U 

A  COMPOSITE  COVARIAICE  VECTOR.  UPPER  TRIAI6LE 

c 

PF  RLO  D 

A  FILTER  COVARIAICE  VECTOR.  UPPER  TRIAIGLE 

c 

T  DP  U  A  SIMDUTIOI  TIME 

c 

XF  RLO  U  A  FILTER  STATE  VECTOR 
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C  IS  RLO  U  A  STSTEH  TRUTH  STATE  VECTOR 

C  /l»RPAR/ 

C  IRXn  II  O  CURREIT  ROI  lOMBER.  1  TO  IRUIS 

C 

C  LOCAL  VARIABLES 
C 

C  PARAMETERS 


c 

IF 

II 

0 

lOMBER  OF  STATES  II  FILTER  MODEL 

c 

IS 

II 

0 

lOMBER  OF  STATES  II  SYSTEM  TROTH  MODEL 

c 

ITC 

n 

0 

SIZE  OF  COMPOSITE  COVARIAICE  TRIAIGLE,  IC(IC4-l)/2 

c 

ITF 

II 

0 

SIZE  OF  FILTER  COVARIAICE  TRIAIGLE.  IF(IF-M)/2 

c 

ITT 

IH 

0 

lOMBER  OF  VARIABLES  II  TRAJECTORY  VECTOR 

c 

TELF 

II 

0 

lOMBER  OF  ELEMEITAL  FILTERS 

c 

PI 

II 

0 

lOMBER  OF  ELEMEITAL  FILTERS 

c 

telf 

II 

0 

lOMBER  OF  ELEMEITAL  FILTERS 

c 

DDMI 

II 

0 

lOMBER  OF  ELEMEITAL  FILTERS  *  1 

c 

TTELF 

II 

0 

lOMBER  OF  ELEMEITAL  FILTERS  1 

c 

PI 

II 

0 

lOMBER  OF  MEASOREMEITS 

c 

ROOTIIES  OSED 

c 

C3 

EMFIL 

IIITPR 

IIITRI  IIVERT  MMATH 

IMPLICIT  CHARACTER  (A-Z) 

C  •«*  PARAMETERS 

IICLODE  ’SIZES.’ 

C  ««•  ARODMEITS 

IITE6ER  MOOT 
C  ***  COMMOI  VARIABLES 
IICLODE  ’XDIITS.’ 

IICLODE  ’DBOQ.’ 
nCLOK  ’DBOGZ.’ 

IICLODE  ’OSRPAR.’ 

IICLODE  ’ICBLK.’ 

IICLODE  ’MMVAR.' 

C  *•*  COMMOI  VARIABLES 
IICLODE  ’COISTS.’ 

IICLODE  ’EVCITL.' 

IICLODE  ’MOIOFF.’ 

IICLODE  ’USRDTI.’ 

IICLODE  ’ROITIM.’ 

C  ***  COMMOI  VARIABLES 
IICLODE  ’MESDO.’ 

IICLODE  ’MESH.’ 

IlCUm  ’ODMAT.’ 

IICLODE  ’MODSTS.’ 

IICLODE  ’MODFIL.’ 

C  •**  LOCAL  VARIABLES 

IITEQER  I.J.K.ELFT 

REAL  XSO(TELF,IS),  ZFO(TELF,RF} . 

R  KFO(TELF,IF), 

»  PCOCTELF.ITC),  PFO(TELF,BTF) , 

»  PEO(TELF,ITF),  TTO<TELF,ITT) , 

t  EFO(TELF,IF).  AFSO(TELF.IF.IS). 

t  CFOCTELF,IF.IF).  DFO(TELF.IF,IF) . 

*  FFO(TELF,IF.IF),  QFO(TELF,IF,IF) 

C 
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C3< 

C4> 

C 


IF  (  HOOT  .EQ.  1  )  THEI 


C 

C 

C 


C 

C 

C 

C 


6 


C 


C 


35 


«I1ITIALIZE  PROBLEM 

IF  (DB6.GE.2)  WRITE  (XOMISCELF) . ' 
t  (''•♦•••••♦♦•••CALL  IIITPR********ELF-  »M4)')ELF 
IF  ((DBO.GT.O}.AMO.(DBa.LE.2))  WRITE 
k  (0.>C*************CALL  I1ITPR********ELP-  »»,I4)*)ELF 
CALL  IIITPR 

ELSEIF  (  MOOT  .EQ.  2  )  THEI 

tlllTIALIZE  ROI 

IF  (T.LE.. 00001)  THEI 
PARCIT-0 
DO  S  J-l.PI 
PARMI(J)-0.0 
COITIIUE 
EIDIF 

IP  (DB0.GE.2)  WRITE  (X0N15(ELF) . > 


A  <»>*******CAU  IIITRI***ELF,IRDI-  » » ,214)  OELF.IRDI 
IF  ((DBG.GT.O)  .AID.  (DBG.LE.2))  WRITE 
t  <0,'(»'*******CAU  IIITRI***ELF.IRrai-  » » ,214) »)ELF,IHDI 


IF  (DB0.GE.2)  WRITE(XDN06(ELF) . 
ft  '(•••••••••CALL  IIITRI*^^ELF,IHDI-  » • ,214) •)ELF,IRDI 

CAU  IIITRI 

DO  35  ELFT-l.TELF 

PROB(ELFT.l)  -  REALCl.O/TELF) 

COITIIUE 


C 


ELSEIF  (  HOOT  .EQ.  3  )  THEI 
C 

C  fCOPY  VECTORS  TO  WORKIIG  VECTORS  FROM  UST  ITERATIOI 

DO  101  I-1,IS 
XS(I)-XSO(ELF,I) 

101  COITIIUE 

DO  102  I-1,IF 
XF(I)-XFO(ELF,I) 

EF(I)-EFO(ELF,I) 

KF(I)-XFO(ELF,I) 

DO  103  J-1,IF 

CF(I,J)«CFO(ELF,I,J) 

DF(I,J)-DFO(ELF,I,J) 

FF(I,J)-FFO(ELF,I,J) 

QF(I.J)-0FO(ELF,I,J) 

103  COITIIUE 

DO  104  J-1,IS 
AFS(I.J)-AFSO(ELF,I,J) 

104  COITIIUE 

102  COITIIUE 


A-4 


DO  105  i-i.nc 

PC(I)-PCO(ELF.I) 

105 

COITIIDE 

DO  106  I-I.ITF 

PF(I)-PFO(ELF.I) 

PE(I)-PEO(ELF.I) 

106 

COITIlOE 

DO  107  I-1,ITT 

YT(I)-TTO(ELF.I) 

lOT 

C 

COITIlOE 

ELSEIF  (  MOOT  .EQ.  4  )  THEI 

C 

tSkVE  VORXIIG  VECTORS  FOR  lEXT  ITERATIOI 

DO  111  I-I.IS 

XSO(ELF,I)-XS(I) 

111 

COITIlOE 

DO  112  I-I.IF 

XFO(ELF,I)-XFa) 

EFa(ELF,I)-EF(I) 

IFO(ELF,I)«KF(I) 

DO  113  J-1,IF 

CFO<ELF,I,J)-CF(I.J) 

DFO(ELF,I,J)-DF(I.J) 

FFO(ELF,I.J)-FF(I.J) 

QFO(ELF.I.J)-QF(I.J) 

113 

COITIlOE 

DO  114  J-I.IS 
*FS0(ELF,I.J)-4FS(I,J) 

114 

COITIlOE 

112 

COITIlOE 

DO  115  I-I.ITC 

PCO(ELF,I)-PC(I) 

115 

COITIlOE 

DO  117  I-1,ITF 

PFO(ELF,I)-PF(I) 

PEO(ELF.I)-PE<I) 

117 

COITIlOE 

DO  119  I-l.ITT 

TTO(ELF,I)-TT(I) 

119 

C 

COITIlOE 

ELSEIF  (  NODT  .EQ.  5  )  THEM 

!  C 

iVRITE  niAE  DATA  TO  OOTFOT  FILE  BLCITOH 

VRITE(XDM02(TTELF)) 

k 

REAL<T), 

((XFHATd,  J)  .I-1,IF) ,  J-1,TELF) , 

k 

XFM, 

k 

PFm. 

1  k 

PFmiD, 

k 

EM. 

ZH, 

PROB, 

k 

(PARESTdXmi,  J)  ,>1,PI) . 

EA, 

EP. 

ft 

. 

TROPAR 
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c 

ELSEIF  (  HOOT  .EQ.  6  )  THEI 

•SAVE  XF  FOR  HMAE  STATE  ESTIMATES  A  SATE  PF  IITO  MATRIX  FORM 
C  FOR  CALCUUTIOI  OF  REAL  ALFA 
C 

CAU  MMATH(  PFMIS.IF.IF.  'M*.  PF.ITF.l,  >M'.  PF.iTF.l  ) 

C 

CALL  MRATH(  XFMIS.IF.l.  >->.  XF.IF.l.  XF.IF.l  ) 

C 

ELSEIF  (  MOOT  .EQ.  7  )  THEI 
C 

C  tPERFORM  MMAE  CALCOLATIOI  REQOIREO  FOR  EACH  ELEMEITAL  FILTER 
C 

DO  12  I-1,IF 

XFMAT(I.ELF)-XF(I) 

12  COITIIUE 

DO  13  I-1,ITF 

PFMAT(I,ELF)-PF(I) 

13  COITIVOE 

C  # 

C  SCALCOUTE  ZRESO  -  ZF()  -  H(.)  *  XF() 

C  tCALCOUTE  (R>  *  A  «R)/2 

CALL  MMATH(  HXFMHS.1M,1.>->.HFILFM.IM.IF.>*MFMIS,IF.1  ) 

C 

CALL  MMATH(  ZRESM.lM.l.  ZSM.lM.l.  HXFMIS,HH,1  ) 

C  #  -  T 

C  tCALCOUTE  ALPRA(.)  -  H(.)  *  P(,)  *  H(.)  *  R() 

CAU  MMATH(  HPF.IM.HF,  HFILFM.IM.HF.  PFMIS.IF.IF) 

C 

CAU  MMATHCHFILFMT.IF.IM.  *«’  .HFILFM.IM.IF.  ’T’  .HFILFM.IM.IF) 
C 

CAU  MMATHCHPHTF.IM.IM.  HPF.IM.IF.  HFILFMT.IF.IM) 
C 

CAU  MMATH(  ALPHAM.IM.IM.  RFM.IM.IM.  HPHTF.IM.IM) 

C 

C  t  -  1 

C  tCALCOUTE  ALPHAMK.)  -  ALPHAMC.) 

C 

CAU  IITERT  (  ALPHAM.  IM.  ALPHAMI  } 

C 

C  t  T 

C  tCALCOUTE  TRAISPOSE  OF  ZRESMTO  -  ZRESMO 

C 

CAU  MMATH(  ZRESMT.I.IM.  ZRESM.lM.l.  ’T'.  ZRESM.lM.l  ) 
C 

C  tCALCOUTE  EXPOIEITIAL  FOR  DEISITT 

C  t  T  -1 

C  t  RO  «  A  «  R() 

C 

CAU  MMATH(ALPHAMIR.IM.1.’-’,ALPRAMI.IR.IN.’*>. ZRESM.lM.l) 

C 

CAU  MMATH(RTAR.1.1.'-'.ZRESMT.1.IM.’«'.ALPHAMIR,IM.1  ) 

C 

RTAR  -  -RTAR/2.0 
RTART(ELF)  -  RTAR 
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C  tCALCDLATE  SQRT  OF  OETERMIIAIT  OF  ALPHAN. 

C  tCALCOLATE  COEFFICZEIT  (SCALE  FACTOR  PRECEDIIG  EXPOIEITIAL) 

C  tCALCOLATE  THE  DEISITY  -  CQEFFICIEIT  •  EXPOIEITIAL 

C 

CALL  MNATHCDET.IN.ni,  .ALPHAN, IM.IN.  ’D'  .ALPHAN, IN, IN) 

C 

EICOEF(ELF}-  1.0/(((2.0*3.141S9236)**(IN/2))*SQRT(DET(1,1))) 
EXPROB (ELF . 1 ) «EXP (RTAR) *EICOEF (ELF) 

C 

ELSEIF  (  NODT  .EQ.  8  )  THEN 
C 

tPERFORN  CALCOLATIOIS  FOR  NDLTIPLE  MODEL  ADAPTIVE  ESTDUTIOI 
OILT  VHEI  THE  LAST  ELENEITAL  FILTER  HAS  BEEI  COMPLETE  FOR 
THIS  PROPAGATIOI/UPDATE  CTCLE 

tCALCOLATE  THE  COIDITIOIAL  PROBABILITIES  FOR  EACH  ELENEITAL  FILTER 
t  1ST  CALCULATE  THE  DEIOMIIATOR  (-  SON  OF  THE  DEISITIES) 

NIIFLAO-0 
DPR0B>0.0 
DO  136  ELFT-1,TELF 

DPROB-DPROB+PROB (ELFT , 1 ) «EXPROB (ELFT . 1 ) 

136  COITIlOE 

f CALCOUTE  THE  COIDITIOIAL  PROBABILITIES 
DO  140  ELFT-1,TELF 

PROB  (ELFT. l)«PROB(ELFT. 1)*EXPR0B (ELFT, 1)/DPRQB 
IF  (PROB(ELFT.l).LT.NIBPROB)  THEI 
NIIFUG-3 

PROB(ELFT,  D-NUPROB 
EIDIF 
140  COITIlOE 

IF  (NIIFUG.GT.I)  THEI 
DPROB-0.0 
DO  146  ELFT-l.TELF 
DPROB>DPROB+PROB (ELFT , 1 ) 

146  COITIlOE 

DO  148  ELFT-1,TELF 
PROB(ELFT,  1)-PR0B(ELFT,  D/DFROB 
148  COITIlOE 

NIIFLAO-0 

EIDIF 

C 

C  tCALCOLATE  THE  BLEIDED  STATE  ESTIMATES  (XFN) 

C 

CALL  MNATH(  XFN.IF.l, ,XFHAT,IF,TELF, .PROB.TELF.l  ) 

C 

C  tCALCOUTE  THE  VARIAICE  OF  STATES  FROM  EACH  ELENEITAL  FILTER 

C  FROM  THE  BLEIDED  OPTIMAL  STATE  ESTIMATES  k  TRAISPOSE 

DO  160  ELFT-l.TELF 
DO  366  I-1,IF 

DXF(I)-XFMAT(I ,ELFT)-XFN(I) 

DXFMATd  ,ELFT)-DXF(I) 

366  COITIlOE 
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CALL  HMATHC  DXFT.l.IF,  DXF.IF.l.  ’T> ,  DXF.iF.l  ) 

C 

C  tCALCDLATE  COVARIAICE  OF  THE  ELENEITAL  FILTERS  art.  THE  BLEROEO 

C  T 

ESTIMATES,  P  -  DXF  *  DXF 

DO  176  I-l.MF 
DO  177  I-l.IF 
PFM(I .K)-DXF(I)*DXFT(1 ,K) 

177  COITIiUE 

176  COHTIlOE 

rrORI  THAT  COVARIAICE  MATRIX  IITO  A  VECTOR 

CALL  MMATH(PFCOL,ITF.l,’V*.PFM.IF,iF,'V».PFM,IF,IF  ) 

SCALC  BLEIDED  COVARIAICE  PFMM  -  PROS  (ELF)  *PFMM  (ELF)  MATRIX 
DO  180  I-I.ITF 

PFMMVT  ( I ) -PROB  (ELFT .  1 )  •  (PFMAT  (I ,  ELFT) -J-PFCOL  (I )  ) 

PFMMV (I) -PFMMV (I)+PFMMVT (I) 

180  COITIIUE 

SSTRIP  OFF  JUST  OIE  COLUMI  OF  ELEMEITAL  FILTER  U-D  COVARIAICES 
C 
C 

C  iCALC  BLEIDED  COVARIAICE  PFMM  -  PROB(ELF)*PFMN(ELF)  MATRIX 

160  COITIIUE 

CALL  MMATH(PFMM,IF,IF.'-*.PFMHV,ITF,1,»H»,PFIDIV,ITF,1  ) 

C 

C  8C0MPUTE  BLEIDED  FILTER  ERROR  EM  IF  REQUIRED: 

C 

CALL  ENFIL 
C 

C  ICALCULATE  THE  ESTIMATES  OF  THE  PARAMETER 
C 

DO  202  J-l.PI 
PAREST(DUNI.J)-0.0 
DO  201  ELFT-l.TELF 

RAREST  (DUMI ,  J)  -PAREST(DOHI ,  J)+PAREST(ELFT ,  J)  •PROB  (ELFT ,  1 ) 

201  COITIIUE 

PARHI(J)-((PARMI(J)*PARCIT)+PAREST(DWa,J))/(PABCIT+l) 

EA ( J) -PAREST (DUMI , J) -TRUPAR( J) 

EP  ( J) - (PAREST (DUMI , J) -TRUP AR ( J) ) **2 
PARCIT-PARCIT+1 

202  COITIIUE 

DO  190  I«1,ITF 
PFMMV(I)-0.0 
190  COITIIUE 

C 

EIDIF 

C 

C  tRETURI  TO  CALLIIG  PROGRAM  UIIT. 

C 

RETURI 

EID 
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C  OPFILE  ~  OUTPOT  PILE  lAMER  t  lUNBERER 
C 

USAGE  —  CALL  QPFILE(XUI.GEIAKE.FRM.RCL.STATS) 


ABSTRACT  ~  OPFILE  ASSI6IS  lUNBERS  AID  lAMES  TO  THE  OUTPUT  FILES: 

XUI02  ,  CITOM  — >  XU1!02(ELFT)  ,  CiTOII(ELFT)  • 

XUI04  ,  D3T0N  — >  XUN04(ELFT)  .  DSTON(ELFT) 

XUI06  .  MSOFE.!!  — >  XUNOS(ELFT)  .  IIPUT(ELFT)  • 

XUI06  .  MDATA  ~>  XUM06(ELFT)  .  MDATA(ELFT) 

XUI07  ,  Scratch  ~>  XU1!07(ELFT) 

XUI08  .  ERRS  — >  XUM08(ELFT)  .  ERRRS(ELFT) 

XUI03  ,  FLIGHT  ~>  ???  Intura  ??? 

XUI09  ,  OLDEHD  -->  ???  future  ??? 

XUIIO  ,  lEVEID  “>  ???  future  77? 

XUI12  ,  EFGAII  “>  777  future  777 
*  WHERE  >ELFT'  IS  THE  ELENEITAL  FILTER  HUMBER, 

C  •  e.t.  0  <  ELFT  <•■  TELF  +  1. 

C  *  WHERE  TELF  IS  THE  HUMBER  OF  ELEMEHTAL  FILTERS 

C  «  AHD  WHERE  THE  FILE  HAMES  BEGIH  WITH  THE  VALUE  OF  ELFT. 

C 
C 

C  LIHITATIOHS  —  HOHE  KHOWH 
C 

C  HEADER  DATE  ~  93/02/22 
C 


c 

ARGUHEHTS 

c 

X0H02 

IH 

0 

UHIT  HUMBER  FOR  FILE 

"OICHTOM" 

c 

XUH03 

IH 

U 

UHIT  HUMBER  FOR  FILE 

"FLIGHT" 

c 

XnH04 

IH 

U 

UHIT  HUMBER  FOR  FILE 

"OIDSTOM" 

c 

XUH05 

IH 

U 

UHIT  HUMBER  FOR  FILE 

"OlIHPUT" 

c 

XUH06 

IH 

0 

UHIT  HUMBER  FOR  FILE 

"OlMDATA" 

c 

XUHOT 

IH 

0 

UHIT  HUMBER  FOR  FILE 

"01"  SCRATCH  FILE 

c 

FOR  ACCUMUUTIHG  PLOT  DATA 

c 

X0H08 

IH 

U 

UHIT  HUMBER  FOR  FILE 

"OlERRRS" 

c 

10109 

IH 

U 

UHIT  HUMBER  FOR  FILE 

"OLEHD" 

c 

XUHIO 

IH 

0 

UHIT  HUMBER  FOR  FILE 

"HDEHD" 

c 

XUH12 

IH 

U 

UHIT  HUMBER  FOR  FILE 

"KFGAIH" 

c 

TELF 

IH 

u 

HUMBER  OF  ELEMEHTAL  FILTERS 

c 

OOMI 

IH 

US 

TELF  +  1 

C 
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C  conioi  VARIABLES  ACCESSED 
C  /DBOG/ 


c 

TTELF 

II 

DS 

EQDAL  TO  DDMI 

c 

ELF 

II 

D 

IDMBER  OF  ELEMEITAL  FILTER  BEIIG  PROCESSED 

c 

ZDM02 

IK) 

D 

DIIT  IDMBER  FOR  FILE  "__CITOM" 

c 

XDM03 

II 

D 

DIIT  IDMBER  FOR  FILE  "FLIGHT" 

c 

XDN04 

IK) 

D 

DIIT  IDMBER  FOR  FILE  "__DSTOM" 

c 

XDM06 

IK) 

D 

OUT  IDMBER  FOR  FILE  "__IBPDT» 

c 

XDM06 

IK) 

D 

DIIT  IDMBER  FOR  FILE  "__MDATA" 

c 

XDNO? 

IK) 

D 

DIIT  IDMBER  FOR  FILE  SCRATCH  FILE 

c 

FOR  ACCOMDUTIK  PLOT  DATA 

c 

XDM08 

IK) 

D 

DIIT  IDMBER  FOR  FILE  "__ERRRS" 

c 

XDM09 

IK) 

D 

DIIT  IDMBER  FOR  FILE  "OLEID" 

c 

XDM!0 

IK) 

D 

DIIT  IDMBER  FOR  FILE  "IDEID" 

c 

XDM!2 

IK) 

D 

DIIT  IDMBER  FOR  FILE  "KFGAII" 

c 

LOCAL  VARIABLES 

c 

XOI 

II 

D 

VALDE  PASSED  TO  OPFILE  FOR  STARTIIG  FILE  IDMBER 

c 

ELFT 

II 

D 

TEMPORARY  VARIABLE  OF  ELEMEITAL  FILTER  IDMBER 

c 

I 

II 

D 

LOOPIIG  VARIABLE 

c 

FRM 

II 

D 

'FORMAT"  OF  THE  FILE 

c 

RCL 

II 

D 

c 

FIAME 

II 

D 

c 

STATS 

II 

D 

c 

GEIAME 

II 

D 

c 

FIOMA 

II 

D 

c 

FIDN!00 

II 

D 

c 

FIDM!0 

II 

D 

c 

FIDM! 

II 

D 

c 

FBDM 

II 

D 

c 

ELF!00 

II 

D 

c 

ELF!0 

II 

D 

c 

ELF! 

II 

D 

c 

C2 

RODTIIES 

DSED 

lOIE 

C 

IMPLICIT  CHARACTER  (A-Z) 

C  ***  PARAMETERS 

IICLDDE  ’SIZES.’ 

IICLDDE  ’XUHITS.’ 

C  ARGDMEITS 
C  ««*  COMMOl  VARIABLES 
CDEB  *0*  COMMOB  VARIABLES 
IICLDDE  ’DBU6.’ 

IICLDDE  ’DBD02.’ 

IICLDDE  ’DSRDYI.’ 

C  ***  LOCAL  VARIABLES 

IITE6ER  ELFT.DDMI.TTELF.I 
CHARACTER*!!  FRM 
CHARACTER*?  FIAME(DDMI) 

CHARACTER*?  STATS 
CHARACTER*5  GEIAME 
CHARACTER*2  FIDMA(DDMI) 

CHARACTER*!  FIDM!00(DDMI)  ,Flim!0(DDMI)  .FIIMI(DDHI) 
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IITE6ER  FNUM(DI]NI).XUI.RCL 
IITEGER«2  ELFIOO.ELFIO.ELFI 


C 

C3> 


TTELF  -  DUMI 

IF  (TELF  .EQ.  1)  TTELF  -  1 
DO  99  ELFT  -  1. TTELF 
FROM  (ELFT)  -ELFT+XUI- 1 
IF  (ELFT  .IE.  TTELF)  THE! 

100  ELF100-FIIM(ELFT)/100 

ELF10-(FIUM(ELFT) - (ELFIOO* 100) ) /lO 
ELF1-FIOI!(ELFT)-(ELF100*100)-(ELF10*10) 

FIIMIOO  (ELFT)  -CHAR(ELF10O+48) 

IF  (FIUM(ELFT)  .LT.  100)  FIUI!100(ELFT)-’0' 
FMDHIO  (ELFT)-CHAR(ELF10-t-48) 

IF  (FIUII(ELFT)  .LT.  10)  FIOM10(ELFT)-*0' 

FIDHl (ELFT)-CHAR(ELFl+48) 

FHUMA (ELFT) -FIUMIO (ELFT) //FIOMl (ELFT) 

IF  (DBG.GE.2)  WRITE  (XUN15 (ELF) .5994) 
fc  FIOMA(ELFT)  .FIOll  100 (ELFT)  ,FHU1110(ELFT) , 

ft  FIDMl (ELFT) .GEHANE 

IF  ((DBG.GT.0).AID.(DBG.LE.2;)  WRITE  (0.5994) 
ft  FWUMA(ELFT) .FIOMl 00 (ELFT) .FIOMIO(ELFT) . 

ft  FIOMl (ELFT) .GEIAME 

ELSE 

FIOMA(ELFT)-'BL’ 

IF  (TTELF  .EQ.  1)  FIDMA(ELFT)-' * 
FIOMIOO(ELFT)-" 

FIOMIO(ELFT)-* > 

FIOMl (ELFT)- ’ ' 

IF  (DBG.GE.2)  WRITE  (XOMIB (ELF) .5995) 
ft  FIDMA(ELFT),FIDM100(ELFT),FIDM10(ELFT). 

ft  FIOMl (ELFT) . GEIAME 

IF  ((DBG.GT.0).AID.(DBG.LE.2))  WRITE  (0.5995) 
ft  FHOMA  (ELFT)  .  FIOMl 00  (ELFT)  ,  FHOMIO  (ELFT)  . 

ft  FIOMl (ELFT) . GEIAME 

EIDIF 

FRAME  (ELFT)-  FIOMA  (ELFT) //GEIAME 
IF  ((FIOM(ELFT) .GT. 100) .AID. (FIOM(ELFT) .LE.200)) 
ft  XDM15(ELFT)  -  FIDM(ELFT) 

IF  ((FIDM(ELFT) .GT.200) .AID. (FIOM(ELFT) .LE.300)) 
ft  XDM02(ELFT)  -  FIOM(ELFT) 

IF  ((FIOM(ELFT) .GT.400) .AID. (FROM (ELFT) .LE.500)) 
ft  XDM04(ELFT)  -  FIDM(ELFT) 

IF  ( (FROM (ELFT) . GT . 500) . AID . (FROM (ELFT) . LE . 600) ) 
ft  XDM05(ELFT)  -  FMDM(ELFT) 

IF  ( (FMDM(ELFT) . GT . 600) . AID. (FROM (ELFT) .LE . 700) ) 
ft  r0M06(ELFT)  -  FROM (ELFT) 

IF  ( (FIDM(ELFT) . GT . 700) . AID. (FMOM(ELFT) . LE . 800) ) 
ft  XDM07(ELFT)  -  FIDM(ELFT) 

IF  ((FIOM(ELFT) .QT.800) .AID. (FIOM(ELFT) .LE.900)) 
ft  XDM08(ELFT)  -  FROM (ELFT) 

IF  ( (FROM (ELFT) . GT . 12) . AID . (FROM (ELFT) . LE . 20) ) 
ft  XDM15(ELFT)  =  FIDM(ELFT) 


All 


IF  ((DBG.6T.0).An).(DBG.LE.2})  WRITE  (0.S993) 
ft  ELFT.FIUH(ELFT),niAIIE(EtFT),FRM.RCL. STATS 

IF  (TTELF  .EQ.  1)  FRAME (ELFT)-GEIAME 
IF  ((XDH  .OE.  500).AHD.(XDH  .LT.600))  THEM 
IF  (DBG.GE.2)  WRITE  (XUM16(ELF)  . ’(>  *0PEI500") ’) 

IF  ((DBG.GT.0).ARD.(DBG.LE.2)) 
ft  WRITE(0.'("0PEI500")’) 

OPEH  (  UIIT  -  FRDMCELFT) . 

ft  FILE  -  FRAME(ELFT). 

ft  FORM  -  FRM. 

ft  STATUS  -  STATS  ) 

ELSEIF  ((XUH  6E.  700) .AMD. (XUM  .LT.800))  THEM 
IF  {DBG.GE.2)  WRITE  (XUM16(ELF)  . '("0PEH700>  *) ’) 

IF  ((DBG.GT.0).AID.(DBG.LE.2)) 
ft  WRITE(0,'("0PEI700")’) 

OPEH  (  UHIT  -  FHUM(ELFT) , 

ft  FORM  -  ’UHFORMATTEO’. 

ft  STATUS  »  'SCRATCH'  ) 

ELSE 

IF  (0BG.GE.2) 

ft  WRITE(XUM15(ELF),'("0PEH???")') 

IF  ((DBG.GT.0).AHD.(DBG.LE.2)) 
ft  WRITE(0,'("0PEH???")') 

OPEH  <  UHIT  >  FHUN(ELFT), 

ft  FILE  -  FHAME(ELFT). 

ft  FORM  -  FRM, 

ft  RECL  -  RCL, 

ft  STATUS  -  STATS  ) 

EHDIF 

C 

IF  ((DBG.GT.0).AHD.(DBG.LE.2))  TKEH 
WRITE(0.5997)  IRUH.ELFT,GEHAME,FHUN(ELFT)  .FHAHE(ELFT)  ,T 
WRITECO ,5990) (XUH02(I) ,1-1 ,DOMI) 

HRITE(0,5991) (XUM04(I) ,1-1 ,DUHI) 

WRITECO , 5987) (XDM05 (I) , I-l ,DUHI) 

WRITE (0,5992) (XUN06(I) ,1-1 ,DUKI) 

WRITECO , 5989) (XUNOTCI) , I-l ,DUHI) 

WRITECO ,5988) (XUH08CI) ,1-1 ,DDKI) 

WRITECO, 5977) CXUM15CI) ,I-1,DDMI) 

EHDIF 

99  COHTIHUE 

5990  FORMAT  ('XUM02  -  '.3(14)) 

5991  FORMAT  C'XDM04  -  ’.3(14)) 

5987  FORMAT  (’XDH05  -  ’,3(14)) 

5992  FORMAT  (’XUM06  -  ’,3(14)) 

5989  FORMAT  (’XDM07  -  ’.3(14)) 

5988  FORMAT  (’XUM08  -  ’.3(14)) 

5977  FORMAT  (’XDM15  -  ’,3(14)) 

5993  FORMAT  (’5993  ’,14,’  FHUH  -  ’,14,’  FRAME  =  ’,A. 

ft  ’  FRM  -  ’, A,’  RCL  -  ’,14,’  STATS  -  ’,A) 

5994  FORMAT  ( ’6994’ , ’4FHDMA(ELFT)-  ’.A,’  FHDMIOO(ELFT)  > 

ft  ,A,’  //FHOMIO(ELFT)  ’,A,’  //FHUMl (ELFT)  ’ ,A, ’GEHAME’ .A) 

6995  FORMAT  ( ’5996’ . ’5FH0MA(ELFT)-  ’,A,’  FHOMIOO(ELFT)  ’ 

ft  ,A,’  //FHUHIOCELFT)  ’,A.’  //FHUMl (ELFT)  ’, A, ’GEHAME’ .A) 

6996  FORMAT  (’6996’ , ’FHDM(’,I4, ’)  -  ’,14) 
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5997  FORMAT  ('5997  IRUl  -  ',14,’ELFT  -  '.14.'  GEMAME  -  '.A. 
k  '  FIAME{'.I4.')  -  ’.A.'  T  -  '.F7.3) 

5998  FORMAT  ('5998'. 'OPFILE  '.A.'  '.14.'  '.A.'  ’.14) 

5999  FORMAT  ('5999', 'OPFILE  IRDI  T  TI  ’ 

k  'TF  TII  TOOT  TOP  TMAX  TIEXT  TSAVE  TMEAS  ’.A) 

6000  FORMAT  (’6000’.’  ELFT  and  TELF  -  ('.14.'  ’.14.’)’) 

6001  FORMAT  (’6001’.’  TTELF  and  0DMI  -  (’.14.*  ’.14.’)’) 


IF  (DB0.GE.2)  WRITE  (X0M15(ELF) . ’ 
k  (”  RETORI  FROM  OPFILE  ••*••••••••♦*••••*••••”)’) 

IF  ((DBG.GT.0).AID.(DBG.LE.2))  WRITE  (0.’ 
k  (  ”  RETORM  FROM  OPFILE  *•**•••••*•••**••••••••*•  ” ) > ) 

IF  (DBG.GE.2) 
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Appendix  B.  Simulation  Plots 


This  appendix  contains  the  plots  of  the  multiple  model  orbit  estimation  problem 
simulations,  including  plots  for  each  failure  type  as  discussed  in  Section  4.1  and  shown 
in  Table  4.1.  For  each  of  these  failures,  the  Figures  contain  the  following  plots  (Figure 
references  are  given  for  the  first  failure  type  only): 

1.  Figure  B.l:  One  plot  of  the  probability  associated  with  each  elemental  filter. 

2.  Figure  B.2  and  Figure  B.3:  Two  figures,  one  for  range  and  one  for  angle,  each  figure 
containing  five  separate  plots.  Each  separate  plot  is  related  to  an  elemental  filter 
and  graphs  the  data  for: 

•  The  mean  residual. 

•  Mean  residual  d:  Monte  Carlo  calculated  standard  deviations. 

•  Zero  mean  ±  filter-predicted  standard  deviations. 

3.  Figures  B.4  to  B.7:  Four  figures,  one  corresponding  to  each  of  the  four  states.  Each 
figure  consists  of  six  separate  plots.  The  plots  correspond  to  the  blended  and  five 
elemental  filters  and  depict  data  for: 

•  The  mean  state  estimation  error. 

•  Mean  estimation  error  ±  Monte  Carlo  calculated  standard  deviations. 

•  Zero  mean  ±  filter-predicted  standard  deviations. 

4.  Figure  B.8:  One  figure  with  two  plots  relating  to  the  MMAE/Bayesian  parameter 
estimation: 

•  One  plot  of  the  mean  estimation  error,  mean  estimation  error  ±  Monte  Carlo 
standard  deviation,  and  zero  ±  filter-predicted  standard  deviations. 

•  The  second  plot  of  the  true  parameter  and  mean  estimated  parameter  values. 
Plots  occur  sequentially  thereafter  in  such  groups  of  eight. 
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B.l  Measurement  Noise  Bias  at  T  >  25 


Table  B.l.  Orbit  Problem  Failure  Case 
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Failure  Type 
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Figure  B.2.  Case  #1  Range  Residuals 
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Figure  B.3.  Case  #1  Theta  Residuals 
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Figure  B.6.  Case  #1  Theta 
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Figure  B.7.  Case  #1  Theta-rate 
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Figure  B,8.  Case  #1  Parameter 
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B.2  Measurement  Noise  Ramp  during  25  <  T  <  35 


Table  B.2.  Orbit  Problem  Failure  Case  #2 
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Failure  Type 
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Figure  B.9,  Case  #2  Probabilities 
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Figure  B.IO.  Case  #2  Range  Residuals 
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Figure  B.13.  Case  #2  Range-rate 
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Figure  B.14.  Case  #2  Theta 
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Figure  B,15.  Case  ^2  Theta-rate 
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Figure  B.16.  Case  #2  Parameter 
B-19 
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Mean  Estimate  of  Variance  Parameter  (RS+ramp.  1  to  7RS,  25<T<35) 


B.S  Measurement  Noise  Bias  during  25  <  T  <35 


Table  6.3.  Orbit  Problem  Failure  Case  #3 


Error 
Case  # 

Failure  Type 

In  Truth  Model 
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Induced 
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Measurement  Noise 

Rs  —  Rso 

Rs  —  Rso  +  6Rsq 

T  <  25,  T  >  35 
25  <  T  <  35 
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Figure  B.17.  Case  #3  Probabilities 
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Figure  B.18.  Case  #3  Range  Residuals 
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Figure  B.19.  Case  #3  Theta  Residuals 
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Figure  B.20.  Case  #3  Range 
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Figure  B.22.  Case  #3  Theta 
B-26 


Mean  Error, H — Sigma, -F— Predicted  (RS+bias  1  to  7RS,  25<T<35) 


Blended  RIter  3  (Baseline) 


eieu-eiei|i 


eiBU-s)aiu 


aieu’eiam 


eieu-eieqi 


Figure  B.23.  Case  #3  Theta-rate 
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B.4  Measurement  Buis  to  .1  for  T  >  10 


Table  B.4.  Orbit  Problem  Failure  Case  #4 
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Figure  B.25.  Case  #4  Probabilities 
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Figure  B.26.  Case  #4  Range  Residuals 
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Figure  B.27.  Case  ^4  Theta  Residuals 
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Figure  B.28.  Case  #4  Range 
B-33 
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Figure  B.29.  Case  #4  Range-rate 
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Figure  B.31.  Case  #4  Theta-rate 
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B.5  Measurement  Ramp  up  to  .1  for  10  <  T  <  20 


Table  B.5.  Orbit  Problem  Failure  Case  #5 
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Figure  B.33.  Case  #5  Probabilities 
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Figure  B.34.  Case  #5  Range  Residuals 
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Figure  B.37.  Case  #5  Range-rate 
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Figure  B.38.  Case  #5  Theta 
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Figure  B.39.  Case  #5  Theta-rate 
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B.6  Measurement  Bias  to  .15  for  T  >  10 


T&ble  B.6.  Orbit  Problem  Failure  Case  #6 
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Figure  B.41.  Case  #6  Probabilities 
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Figure  B.42.  Case  #6  Range  Residuals 
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Figure  B.43.  Case  #6  Theta  Residuals 
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Figure  B.47.  Case  #6  Theta-rate 
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Figure  B.48.  Case  #6  Parameter 
B-55 


Estimate  of  Measurement  Parameterr  (ZS+bias  of  .15,  ®,  T>25) 


B.7  Measurement  Ramp  to  .15  for  10  <  T  <  25 


Table  B.7.  Orbit  Problem  Failure  Case  #7 
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Figure  B.49.  Case  #7  Probabilities 
B-57 
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Figure  B.50.  Case  ^7  Range  Residuals 
B-58 
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Figure  B.53.  Case  #7  Range-rate 
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Figure  B.54.  Case  #7  Theta 
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Figure  B.55,  Case  #7  Theta-rate 
B-63 
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Figure  B.56.  Case  ^7  Parameter 
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13  ABSTV'.if 


Abstract 

Multiple  model  Kalman  Filter  (KF)  techniques  are  used  extensively  for  Multiple  Model  Adaptive  Estimation  (MMAE), 
Multiple  Model  Adaptive  Control  (MMAC),  and  DisUibuted  Kalman  Filter  (DKF)  applications  to  determine  Bayesian-blcnded 
optimal  estimates  of  states,  uncertain  parameters,  and  optimal  conUol  signals.  Multiple  model  methods  are  used  for  sensor 
management.  Failure  Detection  and  Isolation  (FDD.  and  other  Ciuidance  and  Control  (G&C)  applications.  A  simulation  tool 
called  the  Multiple  Model  Simulation  for  Optimal  Filter  Evaluation  (MMSOFE)  has  been  in  this  research.  MMSOFE  is  based 
on  the  well-benchmarked  single  Kalman  filter  tool  called  Multimode  Simulation  for  Optimal  Filter  Evaluation  (MSOFE). 
MMSOFE  is  a  highly  portable  and  versatile  multiple  and  single  Kalman  filter  evaluation  tool.  It  is  capable  of  performing 
simulations  with  one  filter  or  up  to  98  elemental  filters  in  a  multiple  model  adaptive  filter  structure.  It  can  be  adapted  easily 
for  other  multiple  model  applications.  MMSOFE  was  applied  to  failure  detection  and  isolation  of  measurement  jamming-  and 
spoofing-type  failures,  similar  to  jamming  and  spoofing  of  a  Global  Positioning  System  (GPS).  A  satellite  orbit  estimation  test 
case  was  used. 
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